1D modelling of an alpha type Stirling engine

– In this paper, a 1-D model of an half alpha type ‘‘double-acting’’ Stirling engine (made up of two double-acting pistons, two hot heat-exchangers, two cold heat-exchangers and a common recuperator-regenerator for the two separate gas circuits, all that giving two cycles with 180 (cid:2) phasing) is presented. In this architecture, two cylinders, containing, each one, a double-acting 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 ﬂows in the half-engine, under energy transfers. A ﬁnite-volume 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 inﬂuence 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.


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 moderate-temperature parabolic-trough 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 mono-dimensional model, giving the instantaneous speed and pressure of the working gas, derived from the Andersen's one [1].
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 double-acting model, presents a symmetrical configuration around the crankshaft.
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][3][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).

Mono-dimensional model derived from Andersen model
This model allows a mono-dimensional 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)    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).

Mass conservation:
The mass assessment for the 2n numbered main volume is written as follows: 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.
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: For perfect gas, internal energy is a function of the temperature T: U = mAEc v AET.
We thus obtain the expression of the derivative of the temperature in each volume 2n.
The interface temperatures are calculated as follows: Momentum equation: Speeds at the interfaces 2n + 1 are calculated from the equation of the following momentum assessment: 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: with kf is the friction factor and q the gas density.
The pressure at the interface 2n + 1 can be expressed by: 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: The pressures are given here using the differential equation of perfect gases: and

Results of simulation with five main cells
In this first elementary step of simulation, we consider: j five volumes, j constant characteristics of working gas (air): k = 36.10 À3 W/mK; l = 25.10 À6 Pa s, j the exchangers used are with plates: -width of the plates: 200 mm; -heat transfer coefficients are, here, supposed constants [9]: h h = h k = h r = 900 W/m 2 K; -friction coefficient, k f = 0.04 [1,16].
j compression and expansion cylinder volumes are supposed adiabatic, j 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 (p-V) cycles in compression and expansion spaces are plotted only for the last cycle (diagrams (pe-Ve) and (pc-Vc) for the 20th cycle) in Figure 5.
The results obtained with this model, are presented in Table 1.
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): The new results obtained with the Matlab program are presented in Table 2.

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.
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 g i, max are as follows: j For a 20 cells recuperator P i,max = 6625 W and g i,max = 5.4%, j For a 10 cells recuperator P i,max = 6562 W and g 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.

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 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 smooth-plates heat exchangers on the performances of the engine at a speed of 300 rpm.
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.

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 DVw (Figure 9).
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. dQe 2n ¼ À hl 2n A 2n Tw 2n À Te 2n ð Þdt; ð15Þ The heat balance for the element of wall volume DV W allows us to determine the heat stored in the wall: 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: Results of simulation for: -a speed of 300 rpm -an initial pressure of 22.7 bars Table 4. Results of simulation -one-dimensional model, five cells in the recuperator exchangers with smooth plates.

Quantities
Values 5.41    -aluminium walls, for a linear k 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 are presented in Table 5.
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.
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 steady-operating point.
The Table 6 gives the results of simulation under the assumption of the constant wall temperatures in the exchangers.

Study of the combined influence of the out-ofphasing between the two pistons and the initial pressure
The obtained results for an half-engine 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 out-of-phasing 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.

Addition of the thermal and mechanical losses d 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 Dp 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], 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 leak-flow is substracted or added to the working gas flow depending on Dp sign. d 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]: 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. d 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][18][19][20]. First, friction can occur in the Table 5. Results of simulation -one-dimensional model, 10 cells in the recuperator.

Quantities
Values 6.01   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: where, p mf is the average pressure of friction in bar, given by the following equation: 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: The losses by mechanical friction for this alpha type Stirling engine, with double-acting piston, and crankshaft-connecting rod kinematics, are estimated at 2.7 kW.

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 1-D 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 (1-D). 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 outof-phasing of 120°.
This 1-D 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.