Issue 
Int. J. Simul. Multisci. Des. Optim.
Volume 5, 2014



Article Number  A07  
Number of page(s)  11  
DOI  https://doi.org/10.1051/smdo/2013019  
Published online  06 February 2014 
Article
1D modelling of an alpha type Stirling engine
^{1}
EPFEcole d’Ingénieurs, 21 Bd. Berthelot, 34000
Montpellier, France
^{2}
Laboratoire d’Energétique, de Mécanique et d’Electromagnétisme, Université Paris Ouest NanterreLa Défense, 50, rue de Sèvres, 92410
Ville d’Avray, France
^{*} email: nadia_martaj@yahoo.fr
Received:
10
July
2013
Accepted:
13
November
2013
In this paper, a 1D model of an half alpha type “doubleacting” Stirling engine (made up of two doubleacting pistons, two hot heatexchangers, two cold heatexchangers and a common recuperatorregenerator for the two separate gas circuits, all that giving two cycles with 180° phasing) is presented. In this architecture, two cylinders, containing, each one, a doubleacting piston, are combined in order to operate like two parallel Stirling systems in opposition of phase. This model, derived from Andersen’s model, is used to describe the compressible flows in the halfengine, under energy transfers. A finitevolume numerical method is applied for the equations of energy, mass and momentum assessment. This modelling was carried out using the Matlab/Simulink software. The results of this dynamic modelling relate to one half engine. We obtain the evolution of the physical parameters (density, pressure and temperature), and the (p, V) cycle. The influence of the various assumptions was studied. A parametric study was carried out in order to obtain the optimal values of the geometry of the engine and its ideal speed of operation.
Key words: Energetical analysis / Irreversibilities / Optimisation / Stirling engine, Alpha type, Micro solar power plant
© N. Martaj, P. Rochelle, Published by EDP Sciences, 2014
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
In this study we present the results of a Stirling engine simulation, part of a relatively low cost and autonomous micro solar power plant (MiCST) project, based on the use of simple and robust elements to produce electricity for remote sites in emerging countries. Solar heating could be realized by moderatetemperature parabolictrough receivers. The heat storage could be achieved by sensible heat at moderate temperature as T = 300 °C using materials which will be locally available (sand, clay, etc.).
To elaborate a dynamic modelling of the Stirling engine of a micro solar power plant (Figure 1), we built a monodimensional model, giving the instantaneous speed and pressure of the working gas, derived from the Andersen’s one [1].
Figure 1. Schematic of the micro solar power plant. 
The used architecture of the Stirling engine is presented in Figure 2. This alpha type Stirling engine has two cylinders containing, each one, a double effect piston. The two pistons compress the working fluid (air) and receive the mechanical power on their two faces; this doubleacting model, presents a symmetrical configuration around the crankshaft.
Figure 2. Double effect alpha type Stirling engine with recuperator. 
The originality of this engine is that the regenerator or recuperator is a common exchanger for the two parallel circuits in which heat exchanges are done.
This model uses a numerical method of finite volumes. It uses also the energy and the momentum assessment of the compressible flows, and takes into account the internal irreversibilities, like the pressure losses in the exchangers, and includes the dependence of the transfer and friction coefficients to the local Reynolds number [2–4].
For this dynamic modelling of the alpha type Stirling engine, we use the Matlab/Simulink environment.
The engine is mainly made up of at least five volumes (Figure 3), but much more subdivided volumes if necessary:

volume of expansion (e),

volume of the hot exchanger (h),

volume of heat stocking/destocking (recuperator) (r),

volume of the cold exchanger (k),

volume of compression (c).
Figure 3. Subdivision Diagram of the elementary Stirling engine. 
2. Monodimensional model derived from Andersen model
This model allows a monodimensional description of the engine: each volume, such as defined previously, is discretized into one or several main cells. Overlapping these cells, secondary others (with odd numbers) are defined, each one containing the interface of two successive main cells (Figures 3 and 4)
Figure 4. Representative diagram of 3 successive main volumes of the Andersen model with secondary volumes (hatched areas #2n − 1 and #2n + 1). 
In these volumes, we consider the equations of the mass and energy conservation, applied to the main cells and the equation of the momentum in the secondary cells. The physical parameters and sizes can thus be imposed or calculated in each one of these cells and on their walls.
The notations used in the equations are as follows:

n: sequence number of volumes,

2n: order number of the main cells,

2n + 1: order number of the secondary cells.
For example, in the elementary case, where there is no discretization, there are five main volumes 2n (numbers 2, 4, 6, 8, 10), in Figure 3. Volumes 2 and 10 are, respectively, those of expansion and compression. Volumes 4 and 8 are respectively in contact with the hot source and with the cold source. Volume 6 is the recuperator volume.
The equations considered are those of the mass and energy conservation which are applied to the main cells. The equation of the momentum is applied to the secondary cells (which include the junction of the main cells).
2.1. Description of the model
Mass conservation:
The mass assessment for the 2n numbered main volume is written as follows:(1)
The speeds v is counted positives when they are oriented in the order of volumes.
This equation allows, after integration, to obtain the masses in even volumes (main cells).
The masses related to the secondary cells, including the interfaces, are calculated as the average of those of the overlapped volumes.(2)
Energy conservation:
The first law of thermodynamics, for an opened system, applied to the 2n numbered volume, gives the variation of its internal energy U by:(3)
For perfect gas, internal energy is a function of the temperature T: U = m·c _{v} ·T.
We thus obtain the expression of the derivative of the temperature in each volume 2n.(4)
The interface temperatures are calculated as follows:(5)
Momentum equation:
Speeds at the interfaces 2n + 1 are calculated from the equation of the following momentum assessment:(6)where F _{ADi} represents the artificial dissipation force in volume i [1].
The forces due to the pressure losses in the auxiliary cell 2n + 1 are expressed in the following form:(7)where and with kf is the friction factor and ρ the gas density.
The pressure at the interface 2n + 1 can be expressed by:(8)with S _{c,2n} and S _{c,2n+2} as cross sections of the gas flow, in the cells 2n and 2n + 2.
Representative gas speeds in the main cells are:(9)
The pressures are given here using the differential equation of perfect gases:(10)and(11)
2.2. Results of simulation with five main cells
In this first elementary step of simulation, we consider:

five volumes,

constant characteristics of working gas (air): λ = 36.10^{−3} W/mK; μ = 25.10^{−6} Pa s,

the exchangers used are with plates:

compression and expansion cylinder volumes are supposed adiabatic,

each volume is reduced to a central cell (hot and cold exchanger, recuperator, cylinder of compression and cylinder of expansion).
Initial pressure is 22.7 bars (calculated from the initial conditions)
The recuperator temperature is the average of hot and cold temperature sources.
The (pV) cycles in compression and expansion spaces are plotted only for the last cycle (diagrams (peVe) and (pcVc) for the 20th cycle) in Figure 5.
Figure 5. (p, V) diagrams in expansion and compression volumes. 
The results obtained with this model, are presented in Table 1.
Results of simulation – monodimensional model, five volumes, one cell in the recuperator.
The indicated power of the engine, calculated by this model is 5602 W for a velocity about 300 rev/min. This indicated power was calculated by integration on this last cycle. A relative error of 0.54% was calculated, comparing this value to the power obtained by the energy balance. The total heat provided to the cycle is the sum of the heat provided to hot volume Q _{h} and heat Q _{r}, necessary for the compensation of regeneration losses.
The representative speeds of gas in the main cells depend on the local density of the gas which evolves with the temperature. To take into account this dependence, we expressed these speeds by the following equation, which replaces the equation (9):(12)
The new results obtained with the Matlab program are presented in Table 2.
Results of simulation – onedimensional model, five volumes, one cell in the recuperator, dependence of v on the gas density.
2.3. Effect of the subdivision of the recuperator on engine operation
In this section, a monodimensional modelling is presented using several main cells in the regenerator in order to impose a linear variation of temperature in this exchanger. The results obtained by subdividing the recuperator/regenerator in five main cells, using the previous initial data (including the coefficient of transfer of 900 W/m^{2} K, pressure losses 0.04 and the speed about 300 rev/min) are represented in Table 3.
Result of simulation – onedimensional model, five cells in the recuperator.
In order to analyse the influence of the subdivision of the engine volumes on its calculated performance, calculations were made for a recuperator divided initially into 10 then into 20 main cells, maintaining hot and cold volumes subdivided into five main cells, for a speed of 400 rev/min. The results obtained for the maximal indicated power P _{i, max} and the maximal efficiency η _{i, max} are as follows:

For a 20 cells recuperator P _{ i,max} = 6625 W and η _{ i,max} = 5.4%,

For a 10 cells recuperator P _{ i,max} = 6562 W and η _{ i,max} = 5.2%.
We notice that there is not a great difference between the results obtained for these two configurations (10 and 20 cells at the recuperator). The conclusion is that we can limit a reasonable number of subdivisions (10 in this case) to get faster numerical simulations with satisfactory results.
2.4 Addition of friction coefficient correlations on heat exchangers
With the previous initial data, the results are obtained using constant heat transfer and friction coefficients. But physics constraints and reality tell us to use coefficients related to the Reynolds number of the local flow. The correlations of the Nusselt number and the coefficient of friction, used for the smoothplate heat exchangers [10–15], are:

For the laminar flow:

For the turbulent flow:
In these calculations, we considered the following assumptions:

five cells in the recuperator,

one cell for each other volume (hot and cold exchanger, compression cylinder and expansion cylinder).
Table 4 summarizes the influence of the smoothplates heat exchangers on the performances of the engine at a speed of 300 rpm.
Results of simulation – onedimensional model, five cells in the recuperator exchangers with smooth plates.
The figures which follow were obtained using 10 main cells for the recuperator and 5 for each exchanger (22 main cells to represent the engine). The charging pressure is about 10 bars.
We note, in Figures 6–8, the quasi homogeneity of pressure in the circuit, the strong fluctuations of local temperature, which are contradictory to the “isothermal process” assumption, and the speeds of gas synchronism and their homogeneity, except for the sections near the pistons.
Figure 6. Evolution of the pressure according to the angle of crankshaft and to the cell number. 
Figure 7. Evolution of the working gas temperature according to the angle of crankshaft and to the cell number. 
Figure 8. Evolution of the gas speed according to the angle of crankshaft and to the cell number. 
2.5 Addition of the thermal inertia in the solid walls
In the preceding modelling, the wall temperatures of the heat exchangers and of the recuperator were supposed as constants. We will take into account the inertia of these walls which results in a delay in the change of the temperature of wall during transitory operations [6, 7]. We will quantify this delay and the shift of temperature introduced. We will also study its influence on the performances of the engine.
To take into account the thermal inertia of the solid walls, we establish the heat balance of an element of volume ∆Vw (Figure 9).
Figure 9. Exchanger wall cells. 
The Newton’s law allows us to express the heat quantities exchanged by convection between the element of wall volume and the external fluid as well as the element of wall volume and working gas.(15) (16)
The heat balance for the element of wall volume ΔV _{W} allows us to determine the heat stored in the wall:(17)where L _{X} is the length of the exchanger and n _{X} the number of cells per exchanger.
The wall temperature variation is given by:(18)
Results of simulation for:

a speed of 300 rpm

an initial pressure of 22.7 bars

aluminium walls, for a linear λ between 273 W/mK at 300 K and 237 W/mK at 500 K

a mean velocity of coolant of 5 m/s and Nu = 0.023Re ^{0.8} Pr^{0.4}
Results of simulation – onedimensional model, 10 cells in the recuperator.
These results are obtained with five cells for the each exchanger and 10 cells for the recuperator.
We plot in the Figures 10–12, the changes of the internal temperatures of working gas and those of the walls according to the angle of the crankshaft.
Figure 10. Evolution of the gas and wall temperature in the hot exchanger. 
Figure 11. Evolution of the gas and wall temperature in the cold exchanger. 
Figure 12. Evolution of the gas and wall temperature in the recuperator. 
We notice that the temperatures of the walls (in blue) of the hot and cold exchangers and of the recuperator are sinusoidal with very low amplitude (approximately 1 °C), which validates the assumption of the constant local temperature of the walls along the time for a steadyoperating point.
The Table 6 gives the results of simulation under the assumption of the constant wall temperatures in the exchangers.
Results with constant temperatures in the exchanger walls.
2.5 Study of the combined influence of the outofphasing between the two pistons and the initial pressure
The obtained results for an halfengine allow to locate some optimum conditions for operation of this engine.
A parametric study was carried out with three values of out of phasing and two initial pressures:

15 bars and 90°,

10 bars and 110°,

10 bars and 120°.
We obtain a better efficiency (18%) and a power about 8 kW for an outofphasing of 120° (Figure 13). In these conditions, the pressures are not as high in the engine as for the preceding cases, giving a smoother operation.
Figure 13. Evolution of the indicated power and the indicated efficiency versus the engine speeds for different phasing: 90°, 110° and 120° and initial pressure. 
2.6 Addition of the thermal and mechanical losses
● Losses by leakage
The joints are preferably placed at the cold side of the beta type or gamma Stirling engine to avoid the adverse effects due to the temperature; in our case, each piston will be equipped with, at least, one ring. Let us consider a piston surrounded by such a device. If the difference of the pressures is Δp between the two sides of the piston, the instantaneous rate of the mass flow is calculated, such as in a small gap with parallel walls, as follows [5],(19)with L _{d}, the height of the ring.
These losses are estimated at approximately 1500 W for a gap of 1/100 mm and a ring height of 1 cm. The leakflow is substracted or added to the working gas flow depending on ∆p sign.
● Losses by conduction
These losses are due to a heat transfer by conduction on the hot side towards the cold one of the engine [8]:(20)
For a difference between hot and cold sources temperature equal to 200 °C, a length about 0.7 m, a heat transfer coefficient by conduction about 250 W/m.K and a heat transfer area of (0.0126 m^{2}), these losses are estimated at, approximately, 900 W, in absence of insulation between the exchangers.
● mechanical losses
As for the thermal losses, the mechanical losses must also be taken into account to evaluate the power and the total effectiveness of the engine [17–20]. First, friction can occur in the kinematics seals/rings of the piston and the displacer as well as on the shaft. In addition, the seals of the other moving parts can lead to great forces of friction.
The calculation of frictions is carried out using the following expression:(21)where, p _{mf} is the average pressure of friction in bar, given by the following equation:(22)

N: engine speed rev/min

a = 0.5 constant in bar,

b = 0.1 constant in (bar × 1000)/(rev/min).
These values correspond appreciably to half of those used in the reciprocating internal combustion engines. Expressed in term of a friction torque, the expression becomes:(23)
The losses by mechanical friction for this alpha type Stirling engine, with doubleacting piston, and crankshaftconnecting rod kinematics, are estimated at 2.7 kW.
3. Conclusion
In this paper, an original model of double effect Stirling engine was studied. Its originality lies in the exchanges in the regenerator. This is not only a model of a matrix thermal inertia with back and forth exchanges with the same gas, as in the traditional Stirling engines, but an exchanger with two parallel circuits in which the exchanges are done between two gases through a wall with some thermal inertia.
We studied a 1D model, derived from the Andersen model with finite volumes (volumes are subdivided into cells of calculation). We have introduced variable convective heat transfer and flow friction coefficients. This model, using mass, momentum and energy balances of the compressible flows, takes into account the instantaneous internal irreversibilities like the pressure losses in the exchangers and the thermal inertia of the walls as well as the leakages which were taken into account in the most advanced calculus. Using this model, we can plot the speeds, the temperatures, the pressures and the densities of working gas according to time, and/or angle of crankshaft, and space (1D). We obtain the operating characteristics of the engine (powers, efficiencies).
A parametric study showed that a better efficiency (18%) for a power of approximately 8 kW was obtained for an outofphasing of 120°.
This 1D model is simulated under Simulink environment (with multiple cells).
This solar power plant system is under construction and we have not yet the experimental data at our disposal in order to compare with simulation results. Our simulation is used to evaluate the system behavior under several conditions and obtain its power and thermal efficiency and deduce the “best” configuration.
Acknowledgments
The authors gratefully acknowledge Schneider Electric – leader of MiCST project – and all the partners of the consortium: Barriquand Technologies Thermiques, CEAINES, Cedrat Technologies, Defi Systèmes, Exosun, LEME, LEMTA, Mecachrome France, SAED and Stiral. The authors also gratefully acknowledge the ADEME which greatly supports the MiCST project.
References
 Andersen S. March 2006. Numerical Simulation of Cyclic Thermodynamic Processes, PhD thesis, Department of Mechanical Engineering, Technical University of Denmark. (In the text)
 Bonnet S. 2005. Moteurs Thermiques À Apport de Chaleur Externe : Étude d’un Moteur Stirling et d’un Moteur Ericsson. PhD thesis, Université de Pau et des Pays de l’Adour. (In the text)
 Der Minassians A. December 2007. A. Stirling Engines for LowTemperature SolarThermalElectric Power Generation. PhD thesis, University of California at Berkeley.
 Grosu L, Rochelle P. 2009. Application de la méthode de Schmidt avec régénération imparfaite aux 3 types de moteur Stirling. Nouvelles solutions analytiques, congrès SFT, Vol. 2, Vannes, 26–29 mai. p. 895–901. (In the text)
 Homutescu VM, Dumitraşcu G, Horbaniuc B. 2008. Evaluation of the work lost due to leaks through cylinderdisplacer GAP, COFRET’08. Nantes, Juin. (In the text)
 Martaj N. Décembre 2008. Modélisation Energétique et exergétique, simulation et optimisation des moteurs Stirling à faible différence des températures – Confrontation avec l’expérience, Thèse de l’Université Paris 10, Nanterre. (In the text)
 Stouffs P, Bonnet S, Alaphilippe M. 2002. Etude expérimentale des transferts thermiques et des transformations thermodynamiques dans un petit moteur Stirling, in Elsevier, Actes du Congrès SFT’02, Paris. p. 763–768. (In the text)
 Reader Y. 1983. Stirling Engines. E et F.N. Spon, London and New York. (In the text)
 Urieli I, Berchowitz DM. 1984. Stirling Cycle Engine Analysis. Bristol, UK: Adam Hilger. (In the text)
 Tlili I, Timoumi Y, Ben Nasrallah S. 2007. Thermodynamic analysis of Stirling heat engine with regenerative losses and internal irreversibilities. International Journal of Engine Research, 9, 45–56. [CrossRef] (In the text)
 Andersen SK, Carlsen H, Per Grove Thomsen. 2005. Preliminary results from simulations of temperature oscillations in Stirling engine regenerator matrices. Energy, 1–13.
 Kongtragool B, Wongwises S. 2005. Investigation on power output of the gamma configuration low temperature differential Stirling engines. Renewable Energy, 30(3), 465–476. [CrossRef]
 Robson A. 2005. Development of a computer model to simulate a low temperature differential Ringbom Stirling engine, in Thermo and GFD Modelling of Stirling Machines, Proceedings 12th International Stirling Engine Conference, Durham, p. 350–357.
 Ercan Ataera Ö, Karabulut H. 2005. Thermodynamic analysis of the Vtype Stirlingcycle refrigerator. International Journal of Refrigeration, 28, 183–189. [CrossRef]
 Feidt M, Lesaos K, Costea M, Petrescu S. 2002. Optimal allocation of HEX inventory associated with fixed power output or fixed heat transfer rate input. International Journal of Applied Thermodynamics, 5(1), 25–36. (In the text)
 Ibrahim MB, Zhiguo Z, Rong W, Simon TW, Gedeon D. 2002. A 2D CFD Model of Oscillatory Flow With Jets Impinging on a Random Wire Regenerator Matrix. Piscataway NJ, EtatsUnis: IEEE. (In the text)
 Wang JT, Chen J. 2002. Influence of several irreversible losses on the performance of a ferroelectric Stirling refrigerationcycle. Applied Energy, 72, 495–511. [CrossRef] (In the text)
 Senft JR. 1998. Theoretical limits on the performance of Stirling engines. International Journal of Engine Research, 22, 991–1000. [CrossRef]
 Wu F, Chen L, Wu C, Sun F. 1998. Optimum performance of irreversible Stirling engine with imperfect regeneration. Energy Conversion Management, 39(8), 727–732. [CrossRef]
 Organ AJ. 1997. The Regenerator and the Stirling Engine. Wiley, ISBN: 18606580106. (In the text)
Cite this article as: Martaj N & Rochelle P: 1D modelling of an alpha type Stirling engine. Int. J. Simul. Multisci. Des. Optim., 2014, 5, A07.
All Tables
Results of simulation – monodimensional model, five volumes, one cell in the recuperator.
Results of simulation – onedimensional model, five volumes, one cell in the recuperator, dependence of v on the gas density.
Results of simulation – onedimensional model, five cells in the recuperator exchangers with smooth plates.
All Figures
Figure 1. Schematic of the micro solar power plant. 

In the text 
Figure 2. Double effect alpha type Stirling engine with recuperator. 

In the text 
Figure 3. Subdivision Diagram of the elementary Stirling engine. 

In the text 
Figure 4. Representative diagram of 3 successive main volumes of the Andersen model with secondary volumes (hatched areas #2n − 1 and #2n + 1). 

In the text 
Figure 5. (p, V) diagrams in expansion and compression volumes. 

In the text 
Figure 6. Evolution of the pressure according to the angle of crankshaft and to the cell number. 

In the text 
Figure 7. Evolution of the working gas temperature according to the angle of crankshaft and to the cell number. 

In the text 
Figure 8. Evolution of the gas speed according to the angle of crankshaft and to the cell number. 

In the text 
Figure 9. Exchanger wall cells. 

In the text 
Figure 10. Evolution of the gas and wall temperature in the hot exchanger. 

In the text 
Figure 11. Evolution of the gas and wall temperature in the cold exchanger. 

In the text 
Figure 12. Evolution of the gas and wall temperature in the recuperator. 

In the text 
Figure 13. Evolution of the indicated power and the indicated efficiency versus the engine speeds for different phasing: 90°, 110° and 120° and initial pressure. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.