Higher order modeling of a free-piston Stirling engine: analysis and experiment

This paper focuses on higher order modeling and design of the free-piston Stirling engine (FPSE) based on Ant Colony Optimization (ACO). First, the governing thermodynamics and dynamical equations of the engine have been derived. Then, the design parameters of the engine are selected taking into account the finite heat transfer coefficient (resulting in a fifth-order model) and pressure drop (resulting in a sixth-order model) in the dynamical system and the corresponding differential equations are derived in detail. In the following, the mentioned methods and their performance in modeling the FPSE dynamics are investigated. The simulated results show that the effect of the pressure drop on the places of the closed-loop poles of the system is not significant, while the heat transfer coefficient has a considerable effect on the engine dynamics. Accordingly, a fifth-order model along with ACO algorithm is proposed to justify the FPSE behavior. To validate the presented modeling scheme, the prototype engine SUTECH-SR-1 was experimented. It is found that the values of parameters obtained from the proposed design method are close to those of the experiment. Besides, the presented higher order model predicts the engine behavior with an acceptable accuracy through which the validity of the design technique is affirmed.


Introduction
Free-piston Stirling engines (FPSEs) are one of the novel converters for converting solar energy into other types of energy [1,2]. It was first developed by Robert Stirling in the early 1810s in Scotland [3]. Up to now, many researchers have conducted extensive works to design and analyze these types of engines. High efficiency, self-starting, and long operating life can be noted as one of the major advantages [4,5]; however, despite these advantages, designing and setting up these types of engines has been a major challenge for scientists. The dynamic structure of the FPSEs is such that the motion of the pistons of the system must reach the unstable or marginally stable (in linear analysis) and in nonlinear analysis achieve a stable limit cycle [6,7], but unfortunately this condition involves a lot of complexity in the design procedures of Stirling engines. On the other hand, not only parameters such as power, frequency are significant, but also startup condition must be remarked. So far, various practical methods have been documented in literature to analyze and design the FPSEs. The first method is related to the linear analysis so that the behavior of all engine parameters is assumed linearly. The advantage of such systems is their simplicity and reliability simultaneously [8,9]. Linear methods can be discussed and investigated in two directions. The first approach is employed for linear methods in order to analyze and set up FPSE, which is based on the principle that the system is unstable or marginally stable [10]. Truthfully, in these methods, the dynamic analysis of the engine is based on the places of dominant poles of the closed-loop system, and working condition for starting the engine is that the location of closed-loop poles are either on the imaginary axis of the s-plane or on the right side of imaginary axis. The first analysis was carried out to set up and analyze the FPSE employing the wellknown linear method known as the fourth-order model [11,12], which tries to examine the effects of unknown parameters on the engine behavior. The basis of the fourth-order model is devised by Schmidt et al. (also called Schmidt's theory [12], in which the pressure variations in the entire space of the engine compartment are supposed to be, and the heat transfer coefficient in the chamber are also assumed to be infinite. However, the other linear theory presented is known as the fifth-order model, in which the finite heat transfer coefficient is assumed to make it feasible to analyze the behavior of FPSE more precisely [13]. In the case of the second approach, it is used for linear methods to design these engines. Providing a valid design methodology as a means for predicting the critical parameters of the design procedure is crucial although the fact that the FPSEs are passive. Thus far, the most popular method still used is the first method, and there are still few reports regarding the second method. However, in a recent study, Zare and Tavakolpour-Saleh [14] presented an innovative design based on the fourth-order model to design, which has been able to predict optimally the design parameters of the Stirling engine and led to design and construct these engines at a lower cost. But another concern that should be addressed in the design procedure of the desired FPSE is to provide not only an advanced design methodology, but also a more trustworthy method than the fourth-order model, which will undoubtedly lead to important new insights to design with improved accuracy and more reliable to reach as close as possible to the experimental data. For this purpose, in this practical work, the design procedures are conducted with higher order along with Ant Colony Optimization (ACO).
Apart from the first method, the second method relates to nonlinear analysis. In other words, the behavior of one or more system parameters is presumed to be nonlinear. These methods are reliable; however, due to the consideration of nonlinear behavior of the parameters, they involve many complexities such as the phenomenon of limit cycle for the piston [15,16]. Besides the aforementioned works, one of the current challenges in this field is how to propose a method that can be not only practical in the process of precise design, but also reliable for design. This article is an attempt to ignite a method for the mentioned problems.
Today, most of high-tech engineering problems require to optimize at least three parameters which are time-consuming and sometimes expensive. Indeed, it is mainly the advent of the computer in the early 1960s which spurred many advances in the study of complex systems [1]. Hence, with the rapid advancement in computational technology, bioinspired algorithms play a significant role in scientific problems, particularly in engineering fields. In other words, there are numerous practical algorithms, for instance, Genetic Algorithm (GA) [17], Particle Swarm Optimization (PSO) [18], Simulated Annealing (SA) [19,20], etc., and sometimes hybrid algorithms are also performed [21,22], but unfortunately the use of some novel algorithms compared with GA, and PSO are low, especially in FPSE, although the results were satisfactory. In nature, ants accidently search for food (or a source of food), consequently, go back to the colony and itself leave behind itself a chemical substance, namely, pheromone (which can be easily observed in white color after the rain). Other agents, when they find the appropriate route, occasionally, leave the wanderer and follow the new path. ACO has some advantages such as, no central control on ants, quick, neglecting unimportant factors which do not have considerable effect on the system performance, and more importantly, the behavior of each ant is not predictable and some ants opt a longer pass which enable the system to adapt itself to environmental changes [23,24]. These profound investigations made us eager us to employ ACO. It is worth noting that the proposed method in this work has not been reported and ACO has not been performed on such systems.
In this research, a new method for designing the FPSE is presented based on the fifth-order model incorporating ACO technique for the first time. First, two types of fifth and sixth-order models for justifying the FPSEs behaviors are discussed in detail. Comparing the fifth-order model to the sixth-order method reveals that the assumption of instantaneous pressure difference in the expansion and compression spaces (considered in the sixth-order model) have not considerable effect on the engine performance prediction while considering the finite heat transfer is important which is used in the fifth-order model. Therefore, the fifth-order model is proposed for the prototype FPSE. In this work, the mass and spring of pistons are taken into account as important design parameters. Besides, the design process of the engine is based on the locations of the dominant poles or the operating frequency of the system. With this in mind, the amount of unknown design parameters of the engine are obtained by ACO technique and compared to empirical results. Lastly, a test engine is experimented to verify the presented design method and the proposed model.

Free-piston Stirling engines (FPSEs)
An extensive number of works on the design processes of the free-piston Stirling engines (FPSEs) have received widespread attention along with an expanding list of applications lately [3,25]. FPSEs normally consist of six main parts, including the power and displacer pistons, the springs of power and displacer pistons, expansion and compression chambers, and the regenerator as shown in Fig. 1. It should be noted that the regenerator element can be removed from the FPSE components, but the existence of the regenerator will undoubtedly enhance the thermal efficiency (from now on efficiency) of these types of engines [26]. FPSEs are designed according to the Stirling cycle. The ideal Stirling cycle forms a closed-cycle consisting of four linked states including two isothermal and two constant volume regenerative processes as illustrated in Fig. 2. The ideal Stirling cycle efficiency is equal to the Carnot cycle efficiency (Eq. 1). The Carnot cycle consists of two processes including two isothermal and two isentropic processes. To give a meaningful comparison among the thermodynamic cycles, is that the maximum efficiency belongs to the Carnot cycle [27]. According to Fig. 2, the amount of the obtained work in the ideal Stirling cycle is more than the ideal Carnot cycle. As a result, the operation of expansion engines based on the Stirling cycle can be remarked as one of the beneficial aspects of such engines. Figure 2 demonstrates the working principles of the FPSE as follows [4]: (Process 1-2) isothermal compression; heat transfer to the working fluid (at temperature T h ) from the heat source.

Mathematical background investigation
In this section, the governing equations of the FPSE based on the fifth-and sixth-order models are derived and discussed.

Detailed fifth-order model analysis
In this method, the instantaneous pressure in the compression and expansion chambers are assumed to be equal. The dynamic equations governing of the FPSEs based on Newton's second law of motion can be obtained as follows: The instantaneous volumes of expansion and compression spaces can be calculated as: where V h0 , and V co indicate the unavoidable dead volumes. Accordingly, the instantaneous total volume of the working fluid and its derivatives can be stated as: Using Eqs. (8)- (18), the equation of the instantaneous pressure yields: In actual model, the heat transfer between the expansion and compression chambers of the working fluid inside the engine is not infinite and has a constant value. Additionally, using an energetic control volume method, the instantaneous pressure inside the cylinder engine can be calculated. In this analysis, pressure changes as state variables have been added to Schmidt's theory. The internal energy variations using an energetic control volume approach can be expressed as: where Ḣ represents the rate of enthalpy changes, which is considered to be zero in this research. Hence, the instantaneous pressure dynamics can be stated as: Consequently, in order to analyze the system behavior accurately, the engine state variables must first be determined. Therefore, the engine state variables can be written as: To analyze and design an FPSE, based on state variables, the primary step is to convert the instantaneous pressure relationship linearly to the following: As is evident, calculation of constants c 1 -c 5 is crucial to complete Eq. (25). Hence, the unknown constants can be found by linearizing Eq. (19) involves Eqs. (8)- (18) which led to define the unknown constants. With these (20) x 1 = y, interpretations, the closed-loop matrix of the Stirling engine can be easily expressed by: The eigenvalues of Eq. (26) gives five poles as can be expected according to the system order. Therefore, the fifth-order model analysis results in two dominant and three poles as schematically shown in Fig. 3. More importantly, locations of the closed-loop poles on the s-plane carry valuable information regarding the operating frequency of the FPSE system. The imaginary value of the dominant poles denotes the operating frequency of the FPSE system which results in the marginally stability or instability of the close-loop transfer function of the proposed system.

Detailed sixth-order analysis
As discussed previously, instantaneous pressure in the expansion and compression spaces of the engine has never been considered in the fourth-and fifth-order model analysis for the design scheme of the Stirling engine. Hence, in the sixth-order analysis, it is supposed that the instantaneous pressure inside the engine chamber is not equal. According to the mentioned contents, the dynamic equations governing of the engine will be: The total mass inside the engine cylinder of the FPSE is divided into two sections, including the mass of the expansion and compression spaces as: The amount of gas mass inside each cylinder can be obtained easily using ideal gas law as follows: In order to calculate the instantaneous gas pressure in the chambers, first of all, the derivation from the gas mass equation must be calculated: The mass flow rate between the displacer piston and the cylinder wall can be obtained as [28,29]: where c′ is the gap between the displacer piston and the internal FPSE. The instantaneous pressure rate in an expansion chamber can be easily written as follows: Substituting Eqs. (13)-(15) into Eq. (34) and in accordance with Eq. (33), the instantaneous pressure rate in the expansion space can thus be found as: Regarding the calculation of instantaneous pressure rate in the compression space, the derivative of instantaneous volume in compression space can be calculated as: According to Eq. (32), the instantaneous volume rate in the compression space can be expressed as: By substituting Eqs. (14) and (16) into (37), the instantaneous pressure in the compression space can be achieved as: The state variables in this method are: Consequently, in order to investigate the system behavior in space state, the instantaneous pressure rate must be linearized around the equilibrium state. As a result, the linearization of Eq. (38) results in: It is of great importance to find constants a 1 -a 6 and b 1 -b 6 . Hence, we have: Comparing Eqs. (49) and (50), one can easily see that the instantaneous pressure rate in the compression and expansion chambers is equal to the ambient pressure ( P c0 = P h0 = P 0 ). The initial conditions for the state variables are also considered as follows: Considering the initial conditions, the above linearized equations, and instantaneous pressure rate in the expansion and compression spaces [see Eqs. (49) and (50)], the constants a 1 -a 6 and b 1 -b 6 are: With these explanations, the dynamic equations and the instantaneous pressure rate in the expansion and compression spaces can be stated according to the state variables defined as follows: The closed-loop state-space model of the FPSE system can be obtained as follows: The first issue that should be kept in mind is that the eigenvalues of Eq. (75) express closed-loop poles in the dynamic behavior of the FPSE. Indeed, if we could put the dominant pole of the engine's closed-loop on the imaginary axis or to the right hand side of s-plane, the engine will move in motion. According to Eq. (75), the FPSE system includes six poles which schematically show the locations of the closed-loop poles to the working conditions in the marginally stable system (see Fig. 4).

Comparison of the fifth-and sixth-order models
In "Mathematical background investigation", two methods were investigated comprehensively. Now, it is time to decide whether the proposed methods and the effect of parameters such as, the difference in pressure between the expansion and compression spaces, and the limitation of the heat transfer coefficient can have effect on the dynamic performance of the engine or not. For this purpose, this section will compare the methods presented in detail. Therefore, the location of closed-loop poles on the s-plane according to the two methods mentioned above and the empirical information about the SUTRCH-SR-1 engine will be obtained. Furthermore, the comparison of the effect of the two parameters mentioned in the engine performance will be scrutinized. The characteristics of the SUTECH-SR-1 engine are summarized in Table 1.
Regarding the parameters of the SUTECH-SR-1 and the eigenvalues of the state matrix in Eq. (75), the location of the poles of the closed-loop dynamic system for the sixthorder model is obtained in Table 2 and the location of the obtained poles is shown in Fig. 5. Table 2 and Fig. 5 clearly show that the added poles do not have a significant effect on the dynamics of the FPSE, due to the consideration of the effect of the pressure difference between the expansion and the compression chambers, and the two poles added compare with the fourth-order model, which are at the origin and with the changing design parameters, virtually no changes have been made in the location of them. However, according to the specification of the SUTECH-SR-1 engine and the closed-loop system matrix (Eq. 26), the location of the poles according to the fifthorder model is presented in Table 3. It is worth noting that the location of the poles has been obtained by this method is demonstrated evidently in Fig. 6.
The results of Table 3 and Fig. 6 indicate that considering the limited heat transfer coefficient for the engine is a critical parameter and has a definite influence on the dynamics of the engine. By comparing the two methods of the fifth-and sixth-order methods, it can be concluded that the proposed method with the fifth order can be considered as the preferred design method for research. Because the effect of the limited transmission coefficient has a direct impact on the dynamics of the engine, while in the method with the sixthorder model, considering the difference in pressure, there was no effect on the performance of the engine. From now    on, the fifth-order model will be carried out for the rest of the paper. It should be noted that the amount of predicted operating frequency by sixth order is closer to experimental frequency than fifth-order one.

Ant Colony Optimization (ACO)
Ant Colony Optimization (ACO), similar to the available algorithms in literature, inspired from nature, which are the source of the model, has been applied for solving many sophisticated scientific problems so far. ACO was firstly presented by Marco Dorigo in 1992 [30,31] and is one of the most prominent classes of Swarm intelligence. This algorithm uses methods, including edge selection and pheromone update. Ants are unconscious, and have no sound. In addition, the ants cannot see or hear, but they are able to transmit information using a sense of smell (known as stigmergy [32]). When an ant walks, it leaves behind itself a chemical substance called pheromone which evaporates over the time, and the more pheromones, the shorter the pass (perhaps the shortest) can be reached [33]. In other words, choosing a route to achieve the desired goal is not definitive. In order to apply ACO for a problem, memory, artificial obstacles, and living in a discrete environment must be investigated [34,35], respectively. The first one enables researchers to store the route of the moment. The purpose of the second one is to change the details of the problem to attain appropriate results. Finally, the last implementation highlights that ants cannot live without the colony. Detailed graphical information regarding ACO is depicted in Fig. 7.
Currently, the design processes of the FPSE was conducted with the fifth-order method along with ACO. In the proposed method, four parameters have direct effect on the performance of the FPSE. On the other hand, the piston mass and the stiffness attached to it are considered to be unknown. In fact, the main essential of the FPSE is four parameters, such as the power and displacer pistons, as well as the stiffness of the springs. For this aim, the presented paper is an attempt to study this method of using SUTECH-SR-1 engine results to prove the validity of the proposed method. According to the results given in Table 3, the locations of closed-loop poles of the prototype SUTECH-SR-1 engine are achieved. Accordingly, the places of the desired closed-loop poles for the engine are selected based on desired operating frequency (see Table 3). The desired eigenfunction of the FPSE system based on the locations of desired closed-loop poles can be defined as: where 1 indicates the desired eigenfunction. According to the locations of the closed-loop poles of the fifthorder model (see Eq. 26), the desired eigenfunction can be expressed as:  According to the proposed objective function along with ACO algorithm, four unknown design parameters can be calculated. Finally, the proposed procedure of this study is demonstrated as a flowchart (see Fig. 8).

Results and discussion
As discussed earlier, the FPSEs were modeled through three fundamental equations, including Eqs. (3), (4), and (25) which led to the comprehensive state-space model of the FPSE behavior. This section attempts to discuss the fifthorder model along with ACO algorithm to determine four unknown parameters (see Table 4) and, consequently, investigate the performance of the FPSE via the three coupled differential equations. According to the proposed method, the values of unknown parameters are determined with ACO scheme described in "Ant Colony Optimization (ACO)" to achieve the best possible solution to the design problem in the considered ranges of parameters mentioned in Table 4. Figure 9 elucidates not only the convergence of the objective function, but also the minimum acceptable value of 0.00000456 was achieved after 593 iterations.
As a result, the optimal values of design parameters based on the fifth-order method are given in Table 5. The outcomes obtained from the ACO algorithm are compared to GA method. According to Table 5, the results achieved from the ACO not only are in a good agreement with the experimental data, but also are better than those of the GA data, which indicates a sufficient reason to proceed rest of this work with ACO algorithm. Another issue that should be considered is about the time that is taken by an algorithm. With this in mind, the elapsed time of ACO algorithm was 35 min while for GA was about 38 min which shows that the convergence of ACO is faster than the GA. Besides, the obtained design parameters (see Table 5) can be depicted with error bars. Figure 10 shows the error bars of the obtained design parameters and their experimental results. As shown in Fig. 10, the errors between simulated and experimental results are very insignificant.
After obtaining optimum unknown design parameters, FPSEs behavior can be scrutinized concurrently solving dynamic and thermodynamic equations of the engine. As an example, the instantaneous velocity of piston movements of the power displacement can be performed via Simulink/MATLAB. The first point should be noted that is that whether the proposed method and the obtained design parameters are able to achieve a stable motion or not [7]. Thus, first of all, instantaneous velocity and piston displacement have to be investigated. Next, according to the obtained periodic motion, the piston displacement must be studied. Indeed, due to the FPSE behavior, it is vital to achieve a stable periodic motion, and the only strategy is to generate a periodic motion with a set of parameters in the displacer piston [4,7,16]. The simulated instantaneous velocity of the power and displacer pistons according to the fifth-order model is represented in Fig. 11. According to Fig. 11, the maximum speed of power and displacer pistons, respectively, were found as 0.37 and 0.84 m/s, with a better accuracy, through the proposed technique. Moreover, the motions of the piston and power displacements can be obtained according to the design method (see Fig. 12). With regard to Fig. 12, the maximum strokes for power and displacer pistons were obtained as 0.0053 and 0.012 m, respectively. The instantaneous velocity of power and displacer pistons is quite periodic which results in achieving a stable periodic motion.
Another potential issue in the FPSEs must be investigated is the phase portrait of power and displacer pistons [4] as shown in Fig. 13. Furthermore, the stable motion of the FPSE according to the fifth-order model is also affirmed as in fourth model presented in [14]. As results, the fifthorder model achieves a periodic behavior to set up the engine with the obtained design parameters. Now, as a final point, to show the validity of obtained values, Eqs. (13), (14), (17) and also the calculation of pressure variable state were employed to plot the P-V diagram of the prototype engine, namely, SUTECH-SR-1. The area inside the P-V indicates the obtained work per cycle by the engine, which is 0.276 J as depicted in Fig. 14. Another important conclusion must be drawn is about the operating frequency (68.2 rad/s) which results in the output power of 2.97 W.

Experimental study
This section outlines the development of the FPSE based on the obtained design parameters and its specifications are first described. Next, the prototype engine is tested and the validity of the design technique is affirmed via measurement systems.

FPSE specifications
First, the prototype engine namely SUTECH-SR-1 is constructed as shown evidently in Fig. 15. The desired engine is consisted of two power and displacer pistons, while these pistons made of an iron pipe and an iron bar, respectively. The SUTECH-SR-1 engine does not have the regenerator, and the regenerator efficiency is considered to be zero. Besides, the engine cooling part is composed up of a Teflon chamber, which is fed by a pump mounted in the chamber to cool the engine coolant section. Interestingly, the engine structure of the SUTECH-SR-1 is such that the springs attached to piston are easily interchangeable, and the engine can be started up in different working conditions or with any apt spring stiffness. The other thing roughly the engine is that it is easy to assemble and vice versa. The pitch on the head of the engine has made it possible to change the initial space inside the engine by changing the pistons original location, which makes the engine able to adapt itself to different temperatures. Table 6 offers the summary of the Fig. 8 Flowchart of the analytical procedure specification of SUTECH-SR-1. Readers are directed to the cited papers [12] for more details regarding specification of modern FPSEs.

Measurement apparatus
The accelerometer (Model: MMA7361) is used to calculate the power piston displacement. In order to calculate the speed and the displacement of the pistons, the obtained data from the accelerometer through the NI-USB6009 data acquisition are entered into the computer for further analysis. The velocity of the piston of power is achieved by integrating the accelerating the data. Besides, the displacement of piston is obtained by integrating the speed data. In the next subsection, the outcome from the measurement engine is described in detail.

Measurement of work and power
As discussed in "Measurement apparatus", the speed of the power piston can be measured as shown in Fig. 16. As can be inferred from Fig. 16, the maximum speed of power piston is obtained as 0.32 m/s, while in the simulation model, this value was acquired 0.37 m/s. One can obviously see that the extracted maximum values of power piston speed from the experimental and simulation results are in an acceptable agreement. Besides, according to Fig. 17, the maximum power piston displacement achieved is 0.0048 m, while the simulated power piston movement was 0.0053 m. Based on the results obtained from the speed and motion of power piston, the experimental phase portrait plane can be acquired as depicted in Fig. 18. Finally, based on the empirical results, P-V diagram of the engine is plotted. According to Fig. 19, the value of output work obtained is 0.26 J. In addition, according to the experimental power piston displacement, the practical operating frequency can be calculated easily, which is equal to 67.9 rad/s. Consequently, with multiplying the operating frequency and output work, the output power is obtained 2.7 W. Finally, the comparison between simulation and experimental results is demonstrated in Table 7. As depicted in Table 7, there is a good agreement between the simulation and experimental results.

Conclusions
In the present study, the governing thermodynamic and dynamic equations of the engine were first extracted. Then, the fifth-and sixth-order models were described and discussed in detail. Consequently, the effectiveness of the two methods on the FPSE performance using the locations of the closed-loop poles was compared and also evaluated. The results showed that the effect of the pressure drop using the mass flow rate change (sixth order) did not have a considerable effect on the FPSE performance, while it is supposed that the heat transfer coefficient (fifth-order model) has a profound effect on the FPSE performance. Accordingly, the fifth-order model was selected as the design scheme. The aim of the proposed method is to design the essential parameters of the Stirling engine, which plays a fundamental role in the FPSE system, including the power and displacer pistons and springs attached to the pistons. Next, the experimental study was conducted to evaluate the validity of the design procedure as well as the simulation. The validity of the proposed design procedure as well as the corresponding simulation results is a significant issue that must be addressed at this point. Hence, the simulation outcomes were compared to the measured parameters which are clearly in a good agreement. Lastly, the prototype experimental engine (SUTech-SR-1) is developed and experimented. Then, based on the values of designed parameters, the performance of the engine was investigated and compared to the empirical results. The results could predict the operating frequency of the SUTech-SR-1 as 68.2 rad/s with an error less than 0.4%, without causing any practically meaningful error in the design scheme, compared to the experimental results. Moreover, the engine power of 2.97 W was found from the simulation results which 6% is higher than the engine  power compared to the obtained work. To conclude, the results evidently revealed that the effectiveness of the fifth-order model along with ACO algorithm to predict the unknown design parameters. Although our knowledge about the FPSEs has advanced in the past several years, much more advancement remains to be considered which needs more attention in future studies. Additional studies in the field will likely, regarding proposing novel algorithms, not only    result in improved outcomes, but also reduce the computational time.