Experimental and Dynamic Analysis of a Small-Scale Double-Acting Four-Cylinder α -Type Stirling Engine

: This research studies the double-acting four-cylinder α -type Stirling engine. A numerical model is developed by combining the thermodynamic model and dynamic model to study the engine performance. The pressure values of the working zone calculated using the thermodynamic model are taken into the dynamic model to calculate the forces acting on the mechanism. Then, the dynamic model further calculates the displacement, velocity, and acceleration of the mechanism link to provide the pistons’ displacements for the thermodynamic model. The model is also validated using experimental data obtained from testing an engine prototype. Under a heating temperature of 1000 K, cooling temperature of 315 K, charged pressure of 10 bar, and loading torque of 0.33 Nm, the engine is capable of achieving a shaft power of 26.0 W at 754 rpm. In addition, the thermal properties and the transient behavior of the engine can be further simulated using the validated numerical model.


Introduction
The Stirling engine is a device that harvests thermal energy and converts it into mechanical energy. The Stirling engine works based on the temperature difference being applied to the engine. The temperature difference causes the cyclic compression and expansion of working gas, converting thermal energy to mechanical energy through a transmission mechanism. As a Stirling engine is an external combustion engine, it can be adapted with many different heat sources, such as solar energy, geothermal energy, and even wasted heat. Combining solar energy with solar concentrating power systems [1] can reach an efficiency of 32%. There also exists great interest in the harvest of geothermal energy for heating and the cold energy recovery of liquified natural gas for cooling Stirling engines [2]. As an external combustion engine, there is no combustion or explosion during the operation of the Stirling engine. The Stirling engine facilitates a silent operation. Unlike the internal combustion engine, there is no harmful exhaust gas produced by the Stirling engine during its operation.
There are many different configurations of Stirling engines, such as α-type, β-type, and gamma type. The engine studied in this paper has an α-type configuration. A conventional α-type Stirling engine has two cylinders, with each containing one piston. Only one side of each piston is exposed to the working fluid. In our study, the double-acting α-type Stirling engine has both sides of the piston exposed to the working fluid. This doubles the working zone of the engine. The schematics of the four-cylinder double-acting α-type Stirling engine configuration are shown in Figure 1. The engine has four working zones. Each working zone includes the expansion chamber, heater, regenerator, and cooler and is connected to the compression chamber in the consecutive cylinder. An equivalent α-type Stirling engine would have eight cylinders with eight pistons. It is noted that the double-acting α-type configuration is simpler and has a higher power density. Apart from the engine configuration, Stirling engines also have different transmission mechanisms. For the multiple-cylinder double-acting Stirling engines, the transmission mechanism has to be designed while considering the phase angle between the pistons. Hou et al. [3] developed a thermoacoustic double-acting Stirling electrical generator for cold energy recovery. The pistons are directly connected to linear alternators. The developed engine produces 12.4 kW of electric power under low-temperature-differential conditions, in which the cooling and heating temperatures are 110 K and 303 K, respectively. Féniès et al. [4] studied and optimized a double-acting Stirling engine, in which the conventional pistons are replaced with diaphragms connected to springs. The developed engine can reach the power output of 18 W under the charged pressure of 3 bar. Another study by Wu et al. [5] developed a 3 kW double-acting thermoacoustic Stirling electric generator. A maximum electric power of 1.57 kW is produced with 5 MPa pressurized helium at the operating frequency of 86 Hz. Huo et al. [6] proposed a U-shaped transmission mechanism to drive the 4-cylinder double-acting Stirling engine, using crankshafts and gears. Using a numerical simulation, the optimal designs for crankshafts and flywheel aimed at dynamic balance can be found. Campos et al. [7] carried out thermodynamic optimization for a double-acting Stirling engine with the swashplate mechanism. In this study, it is seen that the swept volume and heat transfer area are important parameters for optimization. Carlqvist et al. [8] proposed an in-line configuration for a four-cylinder double-acting Stirling engine, using crankshaft as its driving mechanism. The combustor used is divided into quadrants to heat up the engine manifolds. The engine can reach up to 15 kW of electric power. The wobble yoke transmission mechanism is proposed by Clucas [9] for the development of a four-cylinder double-acting Stirling engine. The engine produces 200 W at 1500 rpm under the charged pressure of 10 bar. In addition, kinematic optimization is carried out by Kojima et al. [10] for a five-cylinder double-acting Stirling engine using the wobble yoke mechanism. Haywood et al. [11] investigated the effects of the number of seals on the performance of a four-cylinder double-acting Stirling engine with the wobble yoke mechanism. His study shows that only a maximum of two seals on the pistons is needed for optimal engine performance. In this study, the wobble yoke mechanism is chosen to drive the four-cylinder double-acting Stirling engine, which gives the pistons a phase angle of 90 • . This transmission mechanism, as shown in Figure 2 converts the translational motion of the piston to the rotational motion of the output shaft.
As the prototyping process and experiments are time-consuming and expensive, it is necessary to predict the engine performance beforehand. This can be done by numerical simulation. The literature [12] has studied different configurations of multiple-cylinder Stirling engines with numerical simulation and found that the configuration with four cylinders connected in series gives the highest power output per cylinder. Cheng et al. [13] coupled the thermodynamic model and the simplified conjugate-gradient method for the optimization of the double-acting Stirling engine. The adopted thermodynamic model is the non-ideal adiabatic thermodynamic model, which considers the heat transfer between the wall and working fluid, as well as the pressure drops within the engine working zone. After optimization, the indicated power output can be increased from 1062.56 to 1659.72 W. Formosa et al. [14] applied an equivalent electrical network model approach for modeling a double-acting Stirling engine with beam springs attached to membranes instead of pistons. The beam springs are modeled as mass-spring-damper systems. The results of the numerical model are compared to experimental results, and the discrepancy in power output is about 30%. Korlu et al. [15] combined the thermodynamic model of both the double-acting Stirling engine and the gas turbine. The gas turbine exhaust gases are used to heat the Stirling engine; therefore, the combined system has a higher efficiency. As the pistons' displacements are affected by pressure variations of two independent working zones, it is important to include the dynamic model. Furthermore, the forces exerted on all four pistons also affect the wobble yoke mechanism. Thus, all four working zones need to be simulated. Cheng et al. [16] simulated the dynamic behavior of the double-acting Stirling engine using the non-ideal adiabatic thermodynamic model and the dynamic model, which consists of force and momentum balance equations derived for a wobble yoke mechanism. The results show that at the charged pressure of 20 atm, the maximum shaft power is 782 W at 247 rpm. García-Canseco [17] simplified the dynamic model for the wobble yoke mechanism into a mass-spring-damper system equation. However, the model still requires the modelling of heat exchangers and pressure drops. Almajri et al. [18] developed a 3D CFD model for predicting the performance of an α-type Stirling engine. Haikarainen et al. [19] simulated the pressure imbalance phenomena in a double-acting Stirling engine by modeling the piston-ring leakage. The simulations show that very small changes in the piston-ring leakage gap can cause pressure instability in the engine working zone. Hachem et al. [20] used the energetic and exergetic models to simulate a double-acting Stirling engine which takes heat from the exhaust hot gas of an internal combustion engine. To date, the existing research still lacks the experimental study of the four-cylinder double-acting Stirling engine with the wobble yoke transmission mechanism. The existing experimental data are also based on the larger scale double-acting Stirling engine. For example, Clucas [9] developed a double-acting Stirling engine with the wobble yoke mechanism that produces 200 W at 1500 rpm under the charged pressure of 10 bar. Pradip [21] analyzed the performance of a WhisperGen MicroCHP engine, which is also a doubleacting Stirling engine with wobble yoke mechanism. The results of the analysis show that the engine gives 1 kW power output and 8 kW thermal power under the charged pressure of 28 bar using nitrogen as the working fluid. Carlqvist and Gothberg [8] developed a four-cylinder double-acting Stirling engine with each bore size of 60 mm. The engine produces 15 kW under the charged pressure of 10 MPa. Since the present paper is focused on a small-scale, mobile, double-acting Stirling engine, it is quite difficult to directly compare our experimental results with the existing reports. Therefore, we need to establish the prototype engine as well as the experimental system so that the experimental data can be collected for the validation of the numerical model. In this study, the non-ideal adiabatic thermodynamic model is coupled with the dynamic model, which is based on force and momentum balance equations derived for the wobble yoke mechanism. The four-cylinder double-acting Stirling engine prototype is also developed. The experiment platform and data acquisition system have been set up to assess the engine performance. The experimental data are then used to validate the simulation model. The simulation model and experimentation will be discussed in detail in the following section.

Numerical Model
This section will explain the developed numerical model for the simulation of the four-cylinder double-acting Stirling engine. The numerical model couples the dynamic model and the thermodynamic model.

Thermodynamic Model
From the dynamic model, the pistons' displacements can be obtained. Then, the pistons' displacements are taken as the input for the thermodynamic model. Based on the pistons' displacements, the volumes of the expansion and compression chambers can be calculated. Then, the pressure variation can be calculated based on the volume variation. In this model, the regenerator is modeled with multiple nodes to capture the temperature gradient phenomenon. While obtaining the pressure variation, the mass change in each sub-chamber can be obtained using the equation below.
In addition, the transport properties of the working fluid, which are the dynamic viscosity and thermal conductivity, are modeled as functions of temperature and pressure based on empirical equations [22]. This is necessary as the values of transport properties vary within the operating range of the Stirling engine. The values of transport properties are needed for the calculation of heat transfer within the engine. The heat transfer coefficients are calculated based on the literature [23,24]. The temperature value of each sub-chamber can be calculated based on the energy conservation equation as below.
The energy equation for working fluid in each sub-chamber is shown below. Additionally, the energy equation for the regenerator mesh is also shown here to describe the heat transfer between the solid matrix and working fluid. Finally, the working zone average pressure of the next time step can be calculated using the ideal gas equation of state. Using helium as the working fluid, the engine has an operating temperature and pressure ranging above helium's critical temperature and critical pressure. When switching to the real gas equation, it is found that the relative difference is within 2%. Therefore, the ideal gas equation of state is used instead of the real gas equation for simplification of the model.
To consider the pressure losses in each chamber and the pressure loss due to sudden expansion and contraction of flow channels, the pressures of the expansion and the compression chambers can be calculated as below. The friction coefficient is obtained from the literature [25,26].
Upon obtaining the pressure values, these values are used to calculate the forces exerted on the pistons' surfaces in the dynamic model.

Dynamic Model
The dynamic model consists of force and momentum balance equations for each link in the wobble yoke mechanism as already shown in Figure 2. Figure 3 shows the free-body diagrams of each link in the wobble yoke mechanism. Based on the free-body diagrams, the force and momentum balance equations of each link can be deduced. From the thermodynamic model, the pressure values of the expansion and compression chambers of all four working zones are taken to calculate the forces acting on the pistons. Therefore, the force balance equation of the piston is shown below.
where F E i represents the force vector acting on point E of the piston, while F c,i and F e,i−1 denote the forces acting on the pistons, which are calculated based on the pressure values of the compression and expansion chambers, respectively. F b represents the force acting on the piston rod from the bounce chamber. In addition, the term m p i g represents the weight of each piston. Then, m p i and a OE i represent the mass and acceleration of the pistons, respectively. Then, the force and momentum balance equations for the piston links can be deduced and shown below. The center of gravity of the piston links is located at the middle of points D and E. Next, for the wobble yoke, considering link C and link D as one rigid link, the center of gravity for the wobble yoke is located at the middle point of point A and point B.
where In addition, I C i C i+2 and I D i D i+2 are the moment of inertia for link C and link D, respectively.
α C i C i+2 and α D i D i+2 are the angular acceleration of link C and link D, respectively. Finally, the force and moment exerted on the tilted plate, which is connected to the main shaft, need to be considered. The force and momentum balance equations are shown below.
where F O and τ O denote the force and torque exerted on the origin, point O. τ f riction and τ loading are the torque due to friction loss and loading, respectively. Here, the value of τ f riction is expressed as a function of the main shaft's rotation speed. The values of the coefficients c 1 and c 2 are tuned to match the experimental data. Solving the abovementioned force and momentum balance equations, the angular accelerations α at each point can be determined. Then, the rotation angle φ and rotation speed of the main shaft ω can be obtained by integrating the angular acceleration of the mainshaft α. Knowing the rotation speed, ω, the piston displacements can be brought into the thermodynamic model for the calculation of the working fluid's thermal properties.

Performance Assessment
Based on the simulation model, the engine performance can be predicted by calculating the indicated power, shaft power, and total efficiency with . W indicated = ω mean 2π P e dV e + P c dV c (19) .
where ω mean is the cyclic mean rotation speed.

Experimental Setup and Methodology
This section discusses the experimental setup and methodology for testing the performance of the prototype four-cylinder double-acting Stirling engine. Figure 4 displays the four-cylinder double-acting α-type Stirling engine prototype developed in this study. Table 1 provides the baseline operating conditions and geometrical parameters for the double-acting Stirling engine prototype.   Figure 5 exhibits the schematic diagram of the experimental setup. The purpose of the experiment is to assess the performance of the engine prototype. In the experiment, the engine is heated using a thermostat with a thermocouple measuring the temperature within the chamber of the electric heater. On the other hand, the engine is cooled with running water with an average temperature of 315 K. The heating and cooling temperatures are measured using another two thermocouples, which are attached to the heater tube surface and the wall of the water jacket. The engine's flywheel is connected to the torque and speed transducers, which measure the rotation speed and torque of the output shaft. The hysteresis brake is used to apply torque on the engine's output shaft. Then, the speed and torque signals are collected by the data logger. Figure 6 shows the photographs of the experimental setup, including the platform for the engine prototype and the data acquisition system. The engine is heated using a FIBROTHALceramic fiber electric heater manufactured by Chung Tai heater enterprise company. The heating temperature is controlled using a thermostat with a PID controller with the accuracy of ±10 • C. The temperature measurements in the experimental setup are carried out using k-type thermocouples, with the measurement range up to 1000 • C, which is suitable for the operating temperature of Stirling engines. The maximum measurement for the torque sensor is 5 Nm, while the maximum engine speed that can be measured by the speed sensor is 5000 rpm. The accuracies for the torque sensor and speed sensor are ±0.02 Nm and ±5 rpm, respectively. The experimental setup uses the hysteresis brake from MAGTROL company, which is used for adjusting the torque exerted on the main shaft. Finally, the signals obtained are acquired using the FLUKE data acquisition system.

Results and Discussions
This section illustrates the results and discussions of this study. It will include the model validation using experimental results, as well as the engine performance curve.

Model Validation
Based on the experimental results, the developed numerical model is validated. Figure 7 shows the comparison between the experimental and simulation results. It can be seen that both sets of results are very close to each other. The operating conditions and geometrical parameters are as in Table 1. At maximum engine speed, the loading torque of 0.2 Nm is applied and gradually increased. As the loading torque increases, the engine speed decreases. Both curves show that there exist different peaks for maximum power output which correspond to different engine speeds. Based on the simulation results, the maximum shaft power is 24.5 W at the engine speed of 780 rpm, under the loading torque of 0.3 Nm.
As for the experimental results, the maximum shaft power of 26.0 W occurs at 754 rpm, under the loading torque of 0.33 Nm. The difference in maximum shaft power output is 5.8%. Therefore, it can be concluded that the numerical model is validated by the experimental results. Regarding the difference between the experimental and numerical results, the difference is mostly due to the calculation of the mechanical loss. To more accurately capture the engine behavior, the equation for frictional loss can be tuned with more coefficient terms.

Baseline Case Simulation Results
Based on the operating conditions and geometrical parameters in Table 1, the engine transient behavior can be simulated using the developed numerical model. Figure 8 shows the variation of shaft power and total efficiency with engine speed. As we can see here, as the engine speed increases, the total efficiency of the engine decreases. This is because more time is needed for heat transfer at a higher engine speed. When the engine is applied to a maximum loading torque of 0.5 Nm, the shaft power is 2.5 W at 48 rpm. A further increase in the loading torque will cause the engine operation to stop. It is also noted that the efficiency of the investigated engine does not exceed 22%. The low efficiency of the engine is significantly due to the large frictional loss in four-cylinder double-acting Stirling engines, compared to single cylinder engines, due to the increased number of pistons. Furthermore, in a small-scale engine, the distance between the hot and cold end is quite short, leading to heat conduction loss. The regenerator can be lengthened to allow for a lower temperature gradient. This also implies that the engine still requires optimization to further improve its performance. Figure 9 conveys the transient behavior of the engine under the baseline operating conditions, with the loading torque of 0.25 Nm, 0.3 Nm, and 0.35 Nm. The final stable operating speeds are 900 rpm, 780 rpm, and 639 rpm, respectively. It can be seen that given the initial engine speed of 1000 rpm, the engine slowly decreases in speed and then reaches a stable engine speed at around 100 s. As the regenerator is modeled using multiple nodes, it takes a longer time to stabilize the regenerator temperature. The response time for the thermodynamic behavior of the engine is longer compared to the dynamic behavior.  The pressure variation of the engine's working zone at the stable regime for the same case is plotted in Figure 10. It can be seen that the pressure variations of each working zone are in a 90 • phase angle with each other. The pressure varies sinusoidally from about 7.3 bar to 11.7 bar. In Figure 10, at the beginning of each engine cycle, each zone represents the pistons' displacements at different main shaft rotation angles, which are 0 • , 90 • , 180 • , and 270 • . The pressure variations of zone 1, zone 2, zone 3, and zone 4 are represented by a solid line, dashed line, dashed-dotted line, and a thicker solid line, respectively. Figures 11 and 12 show the temperature variation of the expansion and compression chambers, respectively, in each working zone. In the expansion chamber, the temperature varies from 850 K to 1041 K. The overcompression phenomenon can be seen here as the expansion chamber temperature becomes higher than the temperature of the boundary condition. For the compression chamber, the temperature varies from 286 K to 346 K. The temperature variations of zone 1, zone 2, zone 3, and zone 4 are represented by a solid line, dashed line, dashed-dotted line, and a thicker dashed line, respectively.   The effects of applying different values of loss coefficient, c 2 , on the engine performance are displayed in Figure 13. In this particular case, the value of c 1 is constant as 1 × 10 −7 kgm 2 /rad. It can be seen that as the value of c 2 decreases, the maximum shaft power output increases, and the performance curve also shifts upward. The higher the value of c 2 , the larger the mechanical loss effects at high engine speed. The values of loss coefficients can be tuned to match the experimental data. The numerical simulation results in Figure 6 are based on c 2 = 2.5 × 10 −7 kgm 2 s/rad 2 .

Parametric Analysis
Parametric analysis of the engine is carried out using the verified model to study the engine under different boundary conditions. Figure 14 shows the effects of charged pressure on the engine power output and total efficiency. In the figure, the performance curves for 10 bar, 12 bar, and 15 bar are shown using a red line, blue line, and green line, respectively. Solid lines represent the engine's shaft power, while the dashed lines represent the total efficiency. Based on the figure, both the shaft power and total efficiency increase as the charged pressure increases. In addition, the engine speed corresponding to the maximum shaft power also increases. When the charged pressure increases to 15 bar, the shaft power is 45.1 W at the engine speed of 940 rpm. This is because when the charged pressure is increased, the force acting on the piston also increases. In addition, the maximum loading torque at 10 bar is 0.5 Nm, while the maximum loading torque at 15 bar is 0.75 Nm. Beyond the maximum loading torque, the engine stops its operation. Without loading torque, the engine can reach an engine speed of 1347 rpm, 1474 rpm, and 1642 rpm at the charged pressure of 10 bar, 12 bar, and 15 bar, respectively. The engine speed increases as the charged pressure increases because there is still enough heat transfer for the engine even though the charged mass increases.  Figure 15 shows the effects of the regenerator's porosity on the engine power output and total efficiency. The figure shows the engine performance curves for the regenerator's porosity of 0.523, 0.691, and 0.768, represented by the red line, blue line, and green line, respectively. The regenerator's porosity of 0.523, 0.691, and 0.768 corresponds to the mesh number of 120, 100, and 50. It is seen that the baseline case regenerator's porosity, which is 0.523, has the highest shaft power and total efficiency. As the porosity increases, the heat transfer surface area of the wire matrix decreases. In addition, when the porosity increases, the dead volume of the engine also increases. Thus, the decrease of heat transfer surface area and the increase of dead volume lead to the decrease of engine shaft power and total efficiency. When the regenerator's porosity is decreased to 0.768, the maximum shaft power is 16.9 W at 683 rpm. Without loading torque, the engine speed can reach 1347 rpm, 1231 rpm, and 1181 rpm for the regenerator's porosity of 0.523, 0.691, and 0.768. On the other hand, the choice of the regenerator's porosity is also limited by the wire mesh sizes available in the industry.

Conclusions
A numerical model that consists of the thermodynamic model and the dynamic model has been developed and validated. Based on the experimental results, the maximum shaft power of 26.0 W occurs at 754 rpm, under the loading torque of 0.33 Nm. The simulated results give the maximum shaft power of 24.5 W at the engine speed of 780 rpm, under the loading torque of 0.30 Nm. The difference in maximum shaft power is 5.8%. Using the validated numerical model, parametric analysis is carried out to study the engine performance under different operating conditions. When the charged pressure increases to 15 bar, the shaft power is 45.1 W at the engine speed of 940 rpm. Furthermore, the maximum loading torque at 10 bar is 0.5 Nm, while the maximum loading torque at 15 bar is 0.75 Nm. When the regenerator's porosity is decreased to 0.768, the maximum shaft power is 16.9 W at 683 rpm. Without the loading torque, the engine speed can reach 1347 rpm, 1231 rpm, and 1181 rpm for the regenerator's porosity of 0.523, 0.691, and 0.768.   Point on wobble yoke mechanism C Point on wobble yoke mechanism c 1 Friction loss coefficient (kg·m 2 /rad)