*Article* **Optimized Piston Motion for an Alpha-Type Stirling Engine**

## **Robin Masser, Abdellah Khodja, Mathias Scheunert, Karsten Schwalbe, Andreas Fischer, Raphael Paul and Karl Heinz Hoffmann \***

Institut für Physik, Technische Universität Chemnitz, 09107 Chemnitz, Germany; robin.masser@physik.tu-chemnitz.de (R.M.); abdellah.khodja@physik.tu-chemnitz.de (A.K.); mathias.scheunert@mailserver.tu-freiberg.de (M.S.); karsten.schwalbe@physik.tu-chemnitz.de (K.S.); andreas.fischer@physik.tu-chemnitz.de (A.F.); raphael.paul@physik.tu-chemnitz.de (R.P.) **\*** Correspondence: hoffmann@physik.tu-chemnitz.de

Received: 26 May 2020; Accepted: 19 June 2020; Published: 23 June 2020

**Abstract:** The Stirling engine is one of the most promising devices for the recovery of waste heat. Its power output can be optimized by several means, in particular by an optimized piston motion. Here, we investigate its potential performance improvements in the presence of dissipative processes. In order to ensure the possibility of a technical implementation and the simplicity of the optimization, we restrict the possible piston movements to a parametrized class of smooth piston motions. In this theoretical study the engine model is based on endoreversible thermodynamics, which allows us to incorporate non-equilibrium heat and mass transfer as well as the friction of the piston motion. The regenerator of the Stirling engine is modeled as ideal. An investigation of the impact of the individual loss mechanisms on the resulting optimized motion is carried out for a wide range of parameter values. We find that an optimization within our restricted piston motion class leads to a power gain of about 50% on average.

**Keywords:** piston motion optimization; endoreversible thermodynamics; stirling engine; irreversibility; power; efficiency; optimization

## **1. Introduction**

In the 1970s, finite-time thermodynamics evolved in Steve Berry's group as an extension to traditional thermodynamics [1]. The aim of this extension was to describe dissipative heat engines operating in finite time or with finite rates as opposed to the reversible description. Finite-time thermodynamics focuses on irreversibilities within the system in question and incorporates them into the analysis. The goal was not so much to capture every little bit of dissipation occurring, the goal was to incorporate the most dominant dissipation contributions in order to get performance features—like efficiencies—much closer to the observed ones than the reversible treatment would give. New concepts [2–5] were developed and then applied not only to heat engines, but also to chemical processes [6]. Already this early work emphasized the importance of process optimization [7]. Later the field widened and different aspects of finite-time thermodynamics were studied such as its usage for power or efficiency optimization [8–10], the influence of different descriptions of irreversibilities [11–15] and the analysis of a broad range of thermodynamic systems [16–19]. Of particular interest are efficiency considerations. From the seminal work of Curzon and Ahlborn [20], to recent publications concerning non-linear irreversible systems [21–24], stochastic fluctuations [25–27], thermoelectric generators [28], biological processes [29] or general realizability domains [30] a multitude of investigations were carried out.

One of the central question in finite-time thermodynamics is: how large is the minimum necessary dissipation to perform a certain process in a specific time? This question can be treated by optimizing

that process—often using control theory—with respect to a certain performance goal, for instance maximizing the power output or minimizing the entropy production. In particular for heat engines a power optimization can be performed by improving the piston motion. Such piston trajectories have been studied for engines with Otto [31–33], Diesel [34–37], and Miller cycles [38] as well as the special paths needed for light-driven engines [39–44].

In this work we will use finite-time thermodynamics methods to study the performance features for Stirling engines. Stirling engines are considered to be good candidates for the use of waste heat, which occurs in many technical applications and which is often dumped into the environment without use. In such situations it is not so much the saved fuel—which comes more or less for free—it is the economical choice of an appropriate engine size (and the connected capital costs) which requires knowledge about performance features and in particular of the power output of the engines. Investigations into such performance features for Stirling engines and their optimization have already been conducted, see for instance [45–50], where especially the piston motion has been considered [51–53].

Here, our goal is to analyze in particular the power output of a Stirling engine in alpha configuration, for which the piston motion is characterized by two independently acting pistons. We determine possible performance improvements for ideal regeneration by varying its piston movement. Finding its performance optima by using the classical approach based on control theory would be a formidable task. Here, we take a simpler route, based on a parameter optimization of the piston motion with especially chosen smooth functions from the AS class of functions introduced below. The advantage of our approach is the simplicity of its technical implementation and a much reduced numerical effort. This turned out to be very favorable in treating the wide range of cases needed in the application which started this investigation: the recovery of waste heat from machine tools. We note that our method will provide lower bounds for the gain achievable by a control theory based optimization. Our investigation will be based on an endoreversible model incorporating the essential losses due to friction of the moving pistons and the resistances in the heat transport in and out of the engine, as well as the impact of the non-vanishing flow resistance in moving the working fluid through the regenerator.

## **2. Piston Motion Optimization**

The standard harmonic piston motion used in Stirling engine modeling is given by

$$V(t) = V\_{\text{dead}} + \Delta V(1 + \sin(2\pi t/t\_0))/2,\tag{1}$$

where *V*dead and Δ*V* are the dead volume and the displacement, respectively, and where *t*<sup>0</sup> is the period of the motion, which is here chosen to be *t*<sup>0</sup> = 1 s. In our analysis we want to capture the effects of two important variations of this standard piston motion. One is a variation in the piston speed as it travels from its minimum displacement to its maximum displacement and back. The other is a variation of the time a working volume is above its average value *V*dead + Δ*V*/2 compared to the time it spends at volumes below its average value. Both variations influence the different loss mechanisms in subtle ways, where the influence of the piston speed on the piston velocity dependent friction loss is most obvious.

Of course the class of piston motions we consider must be periodic, moreover we will consider only motions for which all time derivatives of the piston position exist and are continuous. Making use of the fact that the cross sectional area of the cylinder is constant we directly specify the cylinder volume rather than the piston position. As a result of these considerations we will use the following newly developed piston motion:

$$V(t) = V\_{\text{dead}} + \Delta V f(t/t\_0; \sigma, \delta), \tag{2}$$

where the function *f*(*x*; *σ*, *δ*) and the dimensionless parameters *σ* and *δ* are explained in more detail below. The function *f*(*x*; *σ*, *δ*) is a composition of two functions *f*<sup>1</sup> and *f*<sup>2</sup>

$$f(\mathbf{x}; \sigma, \delta) = f\_1(f\_2(\mathbf{x}; \delta); \sigma), \tag{3}$$

where *f*<sup>1</sup> and *f*<sup>2</sup> are given by

$$f\_1(y; \sigma) = \left(\sin(2\pi y + \sigma \sin(4\pi y)) + 1\right) / 2 \tag{4}$$

and

$$f\_2(\mathbf{x}; \delta) = \mathbf{x} + \delta(1 - \cos(2\pi\mathbf{x})),\tag{5}$$

respectively.

Below we will refer to this class of piston motions as "adjustable sinusoidal" and will label the corresponding items by "AS". Once the piston motion has been optimized by an appropriate choice of parameters we will call it "optimized sinusoidal" and will label the corresponding items "OS". Note that the standard harmonic motion belongs to the class of adjustable sinusoidal motions and can be regained by the choice *σ* = *δ* = 0. We will label items corresponding to this motion "STD". The advantageous feature that the STD case belongs to the AS class of motions allows for an easy comparison and a continuous transition from the STD case to the OS case.

The influence of the parameters *σ* and *δ* on the piston motion is demonstrated in the following figures. Figure 1 shows the changes induced by a variation of *σ*.

**Figure 1.** Piston volume *V* over time *t* for the standard harmonic motion (STD), and the "adjustable sinusoidal" motion (AS) with *δ* = 0 and different values for *σ*. Negative values for *σ* decrease the maximum piston velocity resulting in shorter times spent at the extreme values of the volume, while positive values for *σ* lead to the opposite effect.

As can be seen, an increasing *σ* leads to a faster piston speed (corresponding to steeper slopes of *V*(*t*)) during the compression and expansion phases, while the speed close to the minimum and maximum position is considerably reduced such that the volume stays close to its extreme values for an extended fraction of the overall time period. The figure also shows the effects of negative *σ*. Here, the rest times of the piston at its extreme positions are reduced compared to the standard harmonic motion. This property of the AS motion class is important as it is well known from earlier work on power optimization by motion control [34] that sometimes it is beneficial to let the piston rest for a while at a extreme position.

The effects caused by variations of the parameter *δ* are shown in Figure 2, where the AS motion is displayed for *δ* = −0.05 and *δ* = 0.05. It is apparent that the fraction of time spend above the mean displacement can be extended considerably by setting *δ* to negative values, and can be shortened by setting *δ* to positive values.

**Figure 2.** Piston volume *V* over time *t* for the standard harmonic motion (STD), and the "adjustable sinusoidal" motion (AS) with *σ* = 0 and different values for *δ*. Negative values for *δ* extend the fraction of time with a volume larger than the mean volume. Positive values lead to the opposite effect.

The above discussed features of the AS motion can be observed within certain limits of the parameters *σ* and *δ*, which have been determined by a systematic search of the parameter space. Hence, we will use the following restrictions on these parameters: −0.13 < *σ* < 0.6 and −0.08 < *δ* < 0.08. Beyond these limits the desirable features decline.

While the above discussion described the dynamics of one of the two pistons of the considered Stirling engine, we will now turn to the overall combined dynamics of both pistons. In the subsequent analysis of the potential power improvements through optimized piston dynamics we will use

$$V\_1(t) = V\_{\text{dead}} + \Delta V f(t/t\_0; \sigma\_1, \delta\_1), \tag{6}$$

$$V\_2(t) = V\_{\text{dead}} + \Delta V f(t/t\_0 + \Delta t/t\_0; \sigma\_2, \delta\_2), \tag{7}$$

where we have introduced the additional parameter Δ*t*. This parameter induces a time shift between the two piston motions. In the standard case of harmonic dynamics this parameter is set by the usual phase shift of −*π*/2 which translates into a Δ*t* = −0.25 *t*0. In summary, we will use the five dimensionless parameters *σ*1, *δ*1, *σ*2, *δ*2, and Δ*t*/*t*<sup>0</sup> in our optimization.

## **3. Endoreversible Stirling Engine**

The Stirling engine model used here is based on endoreversible thermodynamics. While a brief description of its concepts is given below, a more extensive description can be found in these two reviews [54,55]. After early work [56–58] the approach became more formalized and is now used in a variety of applications as for instance dealing with chemical reactions [59–63], engines of different kinds [64–67] and in particular with efficiencies of energy transformation devices [68–73].

In the subsequent subsections the endoreversible modeling of the alpha Stirling engine is shown.

## *3.1. Endoreversible Modeling*

The main aspect of endoreversible modeling is the description of a system as well as the processes occurring in it by specifying its subsystems and the reversible or irreversible interactions between those. The subsystems are divided into (in)finite reservoirs with state variables, and engines, which serve for energy conversion. Finite reservoirs are typically described by their energy *Ei* as a function of its extensities, where *i* denotes the *i*-th subsystem. For each extensity *X<sup>α</sup> <sup>i</sup>* a corresponding intensity *<sup>Y</sup><sup>α</sup> <sup>i</sup>* can be calculated by

$$Y\_i^{\mathfrak{a}} = \frac{\partial E\_i(X\_i^{\mathfrak{a}})}{\partial X\_i^{\mathfrak{a}}}.\tag{8}$$

Here, the superscript *α* specifies the corresponding extensity, e. g. the intensity *Y<sup>S</sup>* is the temperature *T*, which is the corresponding intensity to the extensity entropy *S*. The change in energy of the subsystem can thus be expressed as the sum of changes in extensities times the corresponding intensities of the subsystem:

$$\mathbf{d}E\_{\bar{i}} = \sum\_{a} \mathbf{Y}\_{\bar{i}}^{a} \, \mathbf{d}X\_{\bar{i}}^{a} \,. \tag{9}$$

According to this relation, each flux of extensity *J α <sup>i</sup>* carries an accompanying flux of energy

$$I\_i^{\alpha} = \mathcal{Y}\_i^{\alpha} f\_i^{\alpha}. \tag{10}$$

Infinite reservoirs, on the other hand, are only described by their intensities, which stay constant and do not change when extensity or energy is transferred to or from them.

The second type of endoreversible subsystems are engines, which transfer energy from one carrying extensity to another. Unlike reservoirs, they are not intended to store extensities. Hence, extensities and energy are balanced over all incoming and outgoing fluxes:

$$0 = \sum\_{i} f\_{i,k}^{a} \text{ for all } a \text{ and } \tag{11}$$

$$0 = \sum\_{k,\alpha}^{k} I\_{i,k'}^{\alpha} \tag{12}$$

where *k* serves to differentiate the contact points of subsystem *i*, at which the *I α <sup>i</sup>*,*<sup>k</sup>* and *J α <sup>i</sup>*,*<sup>k</sup>* of the various interactions enter the subsystem. Typically, *k* consecutively numbers the contact points of a subsystem, or denotes the linked subsystem. The latter might be favorable when there are only interactions connecting no more than two subsystems—as it is the case in this paper.

The fluxes themselves are defined by either the requirement of equal intensities of two subsystems *Yα <sup>i</sup>* <sup>=</sup> *<sup>Y</sup><sup>α</sup> <sup>j</sup>* or by transport laws for the transferred extensity or energy. While in the first case the fluxes can technically become infinite in order to instantaneously equalize the intensities, in the latter case often phenomenological relationships are used resulting in finite rates. Of course, energy conservation applies to all interactions. In addition, the other extensities must be balanced, with the exception of entropy, since interactions can be irreversible and therefore generate entropy.

Often it is not necessary to describe the energy carrying extensity of an interaction. Instead a power flux is used describing only the rate of transferred energy or work.

#### *3.2. Stirling Engine*

Using the described aspects of endoreversible modeling, we build the Stirling engine model as shown in Figure 3. The subsystems shown as circles are engines representing the regenerator R and the transmission units T1 and T2. The latter are converting the volume work flux of the stroke into some form of power we do not need to specify here. The regenerator is connected to an entropy reservoir SR and a work reservoir WR. It is also connected to the reservoirs representing the gas in the hot cylinder 1 and the cold cylinder 2 of the Stirling engine, respectively. Those in turn are thermally coupled to a hot heat bath H and a cold heat bath C. The remaining reservoirs are work reservoirs collecting the net power WT from the volumetric processes and the frictional power loss WF of the transmission units, as well as volume reservoirs representing the environment E. All reservoirs and interactions are explained in more detail below.

**Figure 3.** Endoreversible model of the Stirling engine with reservoirs (rectangles), endoreversible engines (circles) and reversible (straight lines) and irreversible (wavy lines) interactions. On the left side the hot cylinder 1 is located with its interactions to the hot heat bath H and a transmission unit T1 while on the right side the cold cylinder 2 is displayed with corresponding interactions and the cold heat bath C. Both are connected by the regenerator R in the middle which interacts with an entropy and work reservoir, SR and WR, respectively. Further reservoirs are work reservoirs WT and WF collecting the net power and friction losses, respectively, from the energy converting engines T1 and T2 as well as volume reservoirs E representing the environment.

## *3.3. The Working Fluid*

The working fluid in cylinder 1 and 2 is described by the equation of state for an ideal gas:

$$pV = nRT\_\prime \tag{13}$$

where *p*, *V*, *n* and *T* are the pressure, volume, mole number and temperature, respectively, and *R* is the gas constant. The caloric equation of state is given by

$$
\mathcal{U}I = \pounds\_V nRT,\tag{14}
$$

where *U* is the internal energy and *c*ˆ*<sup>V</sup>* is the dimensionless specific heat capacity at constant volume. Here, this is chosen to be *c*ˆ*<sup>V</sup>* = 5/2 as usual for di-atomic gases.

Endoreversible modelling of reservoirs gets particularly simple, if one uses the extensities to describe its state, which in our case are the entropy *S*, the volume *V* and the mole number *n* of the fluid. We start from the fluid's entropy

$$S = nR\left(\pounds\_V \ln \frac{T}{T\_0} + \ln \frac{V}{V\_0} - \ln \frac{n}{n\_0}\right) + n \frac{S\_0}{n\_0} \tag{15}$$

with *T*0, *V*<sup>0</sup> and *n*<sup>0</sup> being the reference temperature, volume and mole number, respectively, for the reference entropy *S*0(*T*0, *V*0, *n*0). By solving Equation (15) for the temperature *T* we obtain

$$T(S, V, n) = \left(\frac{V\_0 T\_0^{\mathcal{E}\_V}}{n\_0} \frac{n}{V} \exp\left(\frac{S}{nR} - \frac{S\_0}{n\_0 R}\right)\right)^{\frac{1}{\mathcal{E}\_V}}.\tag{16}$$

Now, the internal energy *U* can be expressed in terms of the extensities *S*, *V* and *n*:

$$\mathbb{E}\mathcal{U} = \mathbb{E}\_V n \mathbb{R} \Gamma(\mathbb{S}, V, n) = \mathbb{E}\_V n \mathbb{R} \left( \frac{V\_0 T\_0^{\mathbb{S}\_V}}{n\_0} \frac{n}{V} \exp\left(\frac{\mathbb{S}}{n\mathbb{R}} - \frac{\mathbb{S}\_0}{n\_0 \mathbb{R}}\right) \right)^{\frac{1}{\mathbb{S}\_V}}.\tag{17}$$

This equation is also called the principle equation of state [74]. According to Equation (8) the pressure *p* and the chemical potential *μ* of the fluid can then be derived from the principle equation of state:

$$p(S, V, n) = -\left(\frac{\partial \mathcal{U}}{\partial V}\right)\_{S, n} = \frac{nR}{V}T(S, V, n),\tag{18}$$

$$
\mu(S, V, n) = \left(\frac{\partial \mathcal{U}}{\partial n}\right)\_{S, V} = \left(\mathcal{E}\_V R + R - \frac{S}{n}\right) T(S, V, n). \tag{19}
$$

For a given working fluid databases typically provide the standard molar entropy *Sm*,0(*p*0, *T*0) at reference pressure and reference temperature. In this case, using Equation (13) and *Sm*,0 = *S*0/*n*<sup>0</sup> the temperature can also be expressed as

$$T(S, V, n) = \left(\frac{RT\_0^{1+\varepsilon\_V}}{p\_0} \frac{n}{V} \exp\left(\frac{S}{nR} - \frac{S\_{m,0}}{R}\right)\right)^{\frac{1}{\varepsilon\_V}}.\tag{20}$$

## *3.4. Heat Transfer and Power Losses*

The heat transfer between the hot heat bath H and the gas reservoir 1 as well as between the gas reservoir 2 and the cold heat bath C is assumed to be Newtonian. Thus, the heat flows *I S* 1,H and *I S* 2,C from the hot and cold heat baths are proportional to the temperature differences between the gas reservoirs and the hot and cold heat baths with temperature *T*<sup>H</sup> and *T*C, respectively:

$$d\_{1, \mathcal{H}}^S = \kappa (T\_\mathcal{H} - T\_1)\_\prime \tag{21}$$

$$I\_{2,\mathbb{C}}^{\mathbb{S}} = \kappa (T\_{\mathbb{C}} - T\_2) \, \tag{22}$$

where *κ* is the heat transfer coefficient, which for simplicity is here chosen to be equal for the hot and cold side. The resulting entropy fluxes into or out of the reservoirs are

$$J\_{1, \text{H}}^S = I\_{1, \text{H}}^S / T\_{1'} \tag{23}$$

$$I\_{2,\mathcal{C}}^{\mathcal{S}} = I\_{2,\mathcal{C}}^{\mathcal{S}} / T\_2. \tag{24}$$

The hot and cold heat baths are modeled as infinite reservoirs so that their temperature remains constant at *T*<sup>H</sup> and *T*C, respectively, regardless of the incoming entropy flux.

In order to capture the mechanical friction due to the piston motion, we make the often used assumption of a power loss *P*<sup>f</sup> proportional to the piston velocity squared. Again, since the cross sectional area of the cylinder is constant, this loss can be expressed using the volume change of reservoir *i* = 1, 2:

*<sup>P</sup>*f,i <sup>=</sup> *<sup>β</sup>V*˙ <sup>2</sup> *<sup>i</sup>* , (25)

where *β* is the mechanical friction coefficient [32]. This power loss is transferred to the work reservoir WF for bookkeeping purposes, from which it is then dissipated to heat and dumped into the environment. The resulting power connected to changes of the volume of reservoir *i* is then given by

$$P\_i = p\_i \dot{V}\_i - P\_{\text{f},i} \tag{26}$$

according to the energy balance equations of the transmission units and is flowing to or from the work reservoir WT.

## *3.5. Ideal Regenerator and Mass Transfer*

The regenerator of the Stirling engine is designed to cool and heat the gas flowing back and forth between the cylinders. A real regenerator acts like a short-term heat storage, which during one cycle alternately absorbs energy and then releases it again. Here, we consider an ideal regenerator. We define it in a way such that no irreversibilities occur in the process of regeneration by requiring that the flowing gas leaves or enters the regenerator with the temperature and pressure of the cylinder it flows to or comes from. In addition, the ideal regenerator does not contain any gas itself.

Such an ideal regenerator can be modeled as an engine. The incoming particle flux is passed on directly to the other gas reservoir. The chemical potential and temperature of the incoming and outgoing particle fluxes as well as the associated entropy fluxes are equal to those of the adjacent reservoirs. Hence, these interactions are reversible. We assume a mass transport that is proportional to the difference of the pressures within the two cylinders. Hence, using a mass transfer coefficient *α* the particle flux through the regenerator is given as

$$J\_{1,\mathbb{R}}^{\boldsymbol{n}} = \boldsymbol{a}(p\_2 - p\_1) = -J\_{2,\mathbb{R}}^{\boldsymbol{n}}.\tag{27}$$

The corresponding entropy fluxes entering or leaving the reservoirs 1 and 2 are then given by

$$J\_{1,\mathbb{R}}^S = S\_{m,1} J\_{1,\mathbb{R}'}^n \tag{28}$$

$$J\_{2,\mathbb{R}}^S = S\_{m,2} J\_{2,\mathbb{R}'}^n \tag{29}$$

respectively, where *Sm*,*<sup>i</sup>* = *Si*/*ni* is the molar entropy of the subsystem *i*. We point out that such coupled fluxes can be combined to a multi-extensity flux [75].

Since these two entropy fluxes generally are not equal, a third entropy flux is needed to maintain entropy balance within the ideal regenerator. This third entropy flux *J S* SR,R reversibly flows into the entropy reservoir SR

$$J\_{\rm SR,R}^S = -J\_{1,R}^S - J\_{2,R}^S. \tag{30}$$

We set the temperature *T*<sup>R</sup> = *T*C, as the cold heat bath is considered to be the environment from which entropy or energy can be taken or dumped into at no cost. Likewise, to fulfill energy conservation within the regenerator, an additional energy flux is needed:

$$P\_{\mathcal{R}} = T\_1 f\_{1,\mathcal{R}}^S + \mu\_1 f\_{1,\mathcal{R}}^n + T\_2 f\_{2,\mathcal{R}}^S + \mu\_2 f\_{2,\mathcal{R}}^n + T\_C f\_{\mathcal{SR},\mathcal{R}}^S. \tag{31}$$

This necessary or excess energy for the ideal regeneration process is accounted for in the work reservoir WR and will enter the overall power output. Finally we note, that kinetic energy and mechanical inertia of the gas and mass leakages have been neglected.

## *3.6. The Dynamics*

From the balance equations for the extensities derived above and the transport laws we obtain a coupled system of differential equations to be integrated

$$
\dot{n}\_1 = \mathfrak{a}(p\_2 - p\_1), \tag{3}
\\
\dot{n}\_2 = \mathfrak{a}(p\_1 - p\_2), \tag{32}
$$

$$\mathcal{S}\_1 = \kappa (T\_\text{H} - T\_1) / T\_1 + \mathcal{S}\_{m,1} \dot{n}\_{1\prime} \tag{33} \\ \text{and} \\ \mathcal{S}\_2 = \kappa (T\_\text{C} - T\_2) / T\_2 + \mathcal{S}\_{m,2} \dot{n}\_{2\prime} \tag{33}$$

$$V\_1 = V\_1(t; \upsilon\_1, \delta\_1), \tag{3.1} \\ V\_2 = V\_2(t + \Delta t; \upsilon\_2, \delta\_2), \tag{34}$$

where all the intensities can be expressed in terms of the extensities. For given parameters *σ*1, *δ*1, *σ*2, *δ*<sup>2</sup> and Δ*t*, the above equations are integrated until the system has reached a steady cyclic operation. Resetting the time we obtain the resulting useful work output per cycle

$$\mathcal{W}\_{\text{out}} = \int\_0^{t\_0} P\_\text{R} + P\_1 + P\_2 \,\text{d}t. \tag{35}$$

The overall friction losses can be calculated by

$$P\_{\mathbf{f}} = P\_{\mathbf{f},1} + P\_{\mathbf{f},2} \,. \tag{36}$$

If *W*out is positive, the Stirling engine provides an average power output

$$P = \mathcal{W}\_{\text{out}} / \mathfrak{t}\_0. \tag{37}$$

Otherwise, energy must be supplied to the engine to maintain its operation.

Finally we define the efficiency of the Stirling engine as ratio of the total work output over the integrated heat flux from the hot reservoir for one cycle

$$\eta = \frac{\mathcal{W}\_{\text{out}}}{\mathcal{Q}\_{\text{in}}} = \frac{\mathcal{W}\_{\text{out}}}{\int\_{0}^{t\_{0}} I\_{1,\text{H}}^{\text{S}} \,\text{d}t} = \frac{\mathcal{W}\_{\text{out}}}{\int\_{0}^{t\_{0}} \kappa (T\_{\text{H}} - T\_{1}) \,\text{d}t}. \tag{38}$$

Based on these quantities we now turn to the optimized operation of the Stirling engine.

#### **4. Results**

Our aim was to determine the potential gains in the average power output which could be achieved by an optimized piston motion as compared to the standard motion. To achieve that aim we optimized the parameters *σ*1, *δ*1, *σ*2, *δ*<sup>2</sup> and Δ*t* of the AS motion numerically based on a Nelder–Mead approach [76]. In general we found that the power output was not very sensitive with respect to the motion control parameters in the sense that one got large power variations for small parameter changes.

We presented the results for the piston motion optimized with regard to the average power output (OS) and compared them to those of the standard harmonic piston motion (STD). The results were obtained for model parameters representing a somewhat typical Stirling engine in the few kW range. The temperatures *T*<sup>H</sup> = 400 K and *T*<sup>C</sup> = 300 K reflected the application area of waste heat usage, while the dead volume *V*dead = 1 L and the displacement Δ*V* = 10 L, leading to a typical value of 0.1 for their ratio. The amount of working fluid was *n* = 1 mol, which with the volumes chosen led to a moderate pressure engine. The above engine parameters were kept fixed at their values in the entire results section.

There are three further parameters, for which we investigated their influence on the engine performance in more detail. These were the friction coefficient *β* determining losses due to mechanical friction, the heat transfer coefficient *κ* determining losses due to finite heat conduction, and the mass transfer coefficient *α* which determines the mass flow rate through the regenerator.

We started our analysis by investigating a case in which the power limiting impacts of these three parameters are negligible. In particular the minimum value for the friction coefficient was set to zero, and the values for *κ* and *α* were chosen such that their further increase would no longer lead to a sizable power increase: *<sup>β</sup>*<sup>0</sup> <sup>=</sup> 0 Js/m6 , *<sup>α</sup>*<sup>0</sup> <sup>=</sup> 100 mol/(s bar), *<sup>κ</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup><sup>5</sup> W/K. With these choices the temperature of the working fluid was always very close to the temperature of the connected heat bath and the pressures in the two cylinders were nearly equal. This minimum value for the friction coefficient and the maximum values for *κ* and *α* are referred to as "base values", and the results obtained for these values are referred to as the "base case".

## *4.1. Optimized Piston Motion: The Base Case*

First, we present the results for the base case. The power output for the OS motion is 962 W and for the STD motion 608 W. Thus, the power output for the OS motion turned out to be about 50 % larger than that for the STD motion. The corresponding dynamics of the state variables are shown in Figures 4–8.

In Figure 4 the volumes of the two cylinders are shown as a function of time.

**Figure 4.** Resulting cylinder volumes *V*<sup>1</sup> and *V*<sup>2</sup> over time *t* for the optimized sinusoidal (OS) motion with base case parameters. For comparison the STD motion is plotted with dashed lines.

The optimal volume dynamics showed a number of interesting and surprising features. It is apparent that both volumes varied more like a trapezoid wave than a harmonic wave, leading to a faster transition between the minimal and maximum volumes. Moreover, the extrema were slightly shifted with respect to those of the STD motion. It turned out that the time shift Δ*t*OS = −0.255*t*<sup>0</sup> differed only little from its STD value Δ*t*STD = −0.25*t*0. The values of *δ*1,OS = 0.0217 and *δ*2,OS = 0.036 were positive indicating a preference for smaller volumes for both cylinders, while the relatively large positive values *σ*1,OS = 0.573 and *σ*2,OS = 0.565 reflected the tendency to the square wave behavior. All four parameters differed considerably from the STD case.

The optimized piston motion led to the entropy and mole number dynamics as shown in Figures 5 and 6, respectively. Both figures showed a high degree of similarity, which in part is due to the entropy being roughly proportional to the mole number. Small differences between them were visible around *t* = 0.6 s which indicates that the specific entropy changed there.

**Figure 5.** Entropies *S*<sup>1</sup> and *S*<sup>2</sup> of the hot and cold cylinder, respectively, over time *t* for the OS motion with base case parameters. For comparison those for the STD motion are plotted with dashed lines.

**Figure 6.** Mole numbers *n*<sup>1</sup> and *n*<sup>2</sup> of the hot and cold cylinder, respectively, over time *t* for the OS motion with base case parameters. For comparison those for the STD motion are plotted with dashed lines.

For the further discussion it is helpful to consider the dynamics of the intensities temperature and pressure in the two cylinders. The temperatures of the working fluid in both cylinders were very close to the bath temperatures due to our choice of the base value for the heat conduction. In Figure 7 we show the difference between the cylinder temperature and the corresponding heat bath temperature for both cylinders. One sees that *T*<sup>1</sup> and *T*<sup>2</sup> showed rich dynamics as a consequence of the heat exchange, the volume changes, and the mass transfer through the regenerator. One interesting feature is that during one cycle the cylinder temperature was sometimes largerand sometimes smaller than the corresponding bath temperature. Thus heat entered and returns from each cylinder in relation to its bath. Another feature was the larger variations of the OS temperatures compared to the STD case. The strongly increased power output compared to the STD case required a larger heat supply from the hot bath (and implicitly a larger heat delivery to the cold bath) and therefore larger temperature differences.

**Figure 7.** Temperatures *T*<sup>1</sup> and *T*<sup>2</sup> of the hot and cold cylinder, respectively, over time *t* for the OS and STD motions with base case parameters. Note that the difference to the corresponding heat bath temperature is shown. The temperatures for the OS motion feature a much stronger variation than for the STD case. The rich dynamical structure is due to the interaction of the volume changes in both cylinders.

The pressure dynamics is displayed in Figure 8. Here, the dominating feature was the absence of any noticeable pressure difference between the cylinders during the whole cycle. Again, this feature was caused by our choice of the mass transport coefficient *α*, allowing a sufficiently large mass flow through the regenerator to almost instantaneously equilibrate the pressures of both sides. This common pressure, which was now a global pressure in the whole Stirling engine, showed a much higher peak for the OS motion than for the STD motion. This peak occurred around *t* = 0.6 s, where both volumes reach their minimal values (see Figure 4).

**Figure 8.** Pressures *p*<sup>1</sup> and *p*<sup>2</sup> of the hot and cold cylinder, respectively, over time *t* for the OS motion with base case parameters. For comparison those for the STD motion are plotted with dashed lines. The curves for *p*<sup>1</sup> and *p*<sup>2</sup> lie on top of each other, for better visibility the *p*<sup>1</sup> curves have been moved up by 0.2 bar.

Based on the common pressure *p* = *p*<sup>1</sup> ≈ *p*<sup>2</sup> the infinitesimal total volume work done by both pistons can be expressed as d*W* = *p* d*V*tot with the total volume *V*tot = *V*<sup>1</sup> + *V*2. Then the gain in power output can easily be understood by looking at the *p*-*V*tot-plot, which is displayed in Figure 9.

**Figure 9.** Pressure *p* = *p*<sup>1</sup> ≈ *p*<sup>2</sup> over total volume *V*tot = *V*<sup>1</sup> + *V*<sup>2</sup> for both the OS and STD motion in base case. The OS motion leads to lower volumes and higher pressures resulting in a higher usable work output.

We first noticed that the total volume traversed a larger range for the OS motion than for the STD motion. Especially for small *V*tot this automatically led to much higher pressure values for the given heat bath temperatures. This opened the route to exhaust the work potentially gainable from the considered cycle by reaching into the "corners" of the ideal Stirling cycle with its constant volume branches.

After establishing the base case we then turned to the investigation of the possible power gains by an optimized piston motion in dependence of the loss terms present. We varied the three coefficients one by one to capture the regions of interest for the power output: The minimum value of *β* and the maximum values of *κ* and *α* were chosen to be the base values and the maximum value of *β* and the minimum values of *κ* and *α* were set by a vanishing power output. When varying one of the parameters, the base values were used for the other parameters.

#### *4.2. Optimized Piston Motion: Friction*

First, we looked at achievable power gains for different friction losses. To that end we varied the friction coefficient *<sup>β</sup>* between zero and about <sup>7</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> Js/m6 , for which all the power produced was dissipated by the friction losses.

Indeed the optimized piston motion led to an increase in the average power output compared to the standard motion, as shown in Figure 10. On average the gain in power output was about 50 %, even though the piston motion was only optimized within the AS class considered. While for small friction the power gain was maximal it declined towards larger values of the friction coefficient. Astonishingly it kept a constant power gain from about *<sup>β</sup>* <sup>=</sup> <sup>4</sup> <sup>×</sup> 105 Js/m6 down to vanishing power, which means that there was strong increase in relative performance.

**Figure 10.** Average power output *P* over varied friction coefficient *β* for both OS and STD motion. The friction coefficient has been increased from zero until no positive average power output was reached.

In Figure 11 we show the efficiencies corresponding to the power optimized motion and the standard motion.

**Figure 11.** Efficiency *η* over varied friction coefficient *β* for both OS and STD motion.

It is interesting to note that for vanishing friction losses the efficiencies both started at about 0.25 which corresponded to the Carnot efficiency for the given heat bath temperatures. Moreover the OS efficiency was initially lower than the STD efficiency while the power output was higher by more than 50 %. This means that the large increase in power needed disproportionately more heat input than the STD case, which for a waste heat application was not so crucial.

#### *4.3. Optimized Piston Motion: Heat Conduction*

The second dissipative process present in the Stirling engine is the finite heat conduction. While for high heat conduction the temperatures in the working volumes were close to the temperatures of the associated heat baths, this was different for a low heat conduction. The heat transfer coefficient *κ* was varied between its base value and nearly zero.

As can be seen in Figure 12, the average power output at high heat conduction for the OS motion reached values of about 150 % of that of the STD motion. With *κ* getting smaller the OS power output decreased faster than the STD power output. When *κ* reached values about 4 kW/K the decay in OS power became stronger, but it stayed always above the STD one. For even smaller heat transfer coefficients the power output decreased towards zero. A similar behavior can be observed for the standard motion, however with a considerably smaller power output. While for the OS motion the average power output became negative for *κ* close to zero, this happened for the standard motion already at *κ* values around 200 W/K.

**Figure 12.** Average power output *P* over varied heat transfer coefficient *κ* for both OS and STD motion. The OS motion reaches average power output values of around 150 % compared to the STD motion as well as positive values of *P* with *κ* close to zero.

The efficiency is shown in Figure 13 as a function of the heat transfer coefficient. It is apparent that the efficiency was mostly higher for the standard motion than for the OS motion. However at smaller values of *κ* it dropped below that for the OS motion. It stayed close to the Carnot value 0.25 over most of the considered *κ* range and dropped off around the values where the power output also decreased. This led to a crossing point of the STD and OS graphs in Figure 13 which was not present in Figure 12. There the STD graph was always below the OS graph and thus the decline of the power output for the STD case to zero at larger values of *κ* did not necessitate such a crossing.

**Figure 13.** Efficiency *η* over varied heat transfer coefficient *κ* for both OS and STD motion. For higher values of *κ* the efficiency of the OS motion is slightly lower than that of the STD motion. For lower *κ* the OS motion leads to better efficiency values.

## *4.4. Optimized Piston Motion: Mass Transport*

The third non-equilibrium process in our model is the mass transport between working volumes 1 and 2. This transport was assumed to be proportional to the pressure difference, with the proportionality constant being the mass transport coefficient *α*. However, even though it was a non-equilibrium transport, it described a reversible exchange of working fluid between the working volumes 1 and 2 through the regenerator. This is due to the regenerator being modeled as ideal which corresponded to a fully reversible operation.

In Figure 14 the average power output is shown as a function of the mass transport coefficient *α* for the range between zero and 10 mol/(s bar). For large values the power output saturated, while for smaller values the power output slowly decreased. Around 1.5 mol/(s bar) it started to fall off towards zero. Initially for large enough values of *α* the gain in the power output of the OS motion compared to that of the standard motion was more than 50 %. The gain increased towards smaller *α* and became infinite in the range where the OS power output was positive while the STD power output was negative.

**Figure 14.** Average power output *P* over varied mass transfer coefficient *α* for both OS and STD motion. The OS motion leads to an increase in *P* of more than 50 % and to lower feasible values of *α* compared to the STD motion.

In Figure 15 the efficiency is shown as a function of the mass transport coefficient *α*. This figure shows efficiencies very close to the Carnot efficiency *η*<sup>C</sup> = 0.25, which is due to the reversible modeling of the regenerator. For very small *α* the finite—but large—heat transfer coefficient *κ* led to negative power output, for which we set the efficiency to zero. The surprising feature here was the steep decline of efficiency once it left the Carnot value level.

**Figure 15.** Efficiency *η* over varied mass transfer coefficient *α* for both OS and STD motion. Despite the equally high efficiency values above *α* > 0.8 mol/(s bar), the OS motion maintains such high values for much lower mass transfer coefficients than the STD motion.

## **5. Conclusions**

In this paper we investigated possible performance improvements for an endoreversible Stirling engine through an optimized piston motion. The motion of the two pistons of the Stirling engine was limited to the smooth adjustable sinusoidal (AS) motion which is controlled by five parameters. The underlying endoreversible model that was build for the Stirling engine takes into account friction losses, irreversible heat transfers as well as the impact of a finite gas flow through the regenerator. The regenerator itself was assumed to be ideal.

For comparison a "base case" was defined without friction losses and where heat and mass transfer coefficients have been chosen such that their further increase would no longer lead to a sizable increase in power output. Already in this case, the piston motion optimization shows that the average power output can be increased by about 50 %. These surprising gains could be achieved by the fact that the optimized piston motion resulted in higher temporal pressure variations within the system. This brought the process closer to the ideal Stirling cycle as is well illustrated in the pressure-volume diagram.

In all cases considered, there was an astonishing gain in average power output. This effect is especially pronounced in cases with "unfavorable" coefficients, that means high friction or low heat and mass transfer coefficients. This extends to conditions, where the Stirling engine with standard piston motion can no longer be profitably operated, but where it can still be operated with good performance using an optimized sinusoidal (OS) piston motion.

The efficiency of the Stirling engine with optimized piston motion was larger or smaller, depending on each case, compared to the standard motion. This is a consequence of the optimization being done with the objective to maximize the average power rather than the efficiency.

The fact that we have used an ideal regenerator has a large impact on the system behavior as well as the optimization of the piston motion. This work was carried out to show the potential gains in power output that can be achieved by optimizing the piston motion. While the results indicate a considerable performance increase for the investigated model with an ideal regenerator, real engines will have non-ideal regeneration. Our future work will therefore focus on a modified endoreversible Stirling engine with a non-ideal regenerator.

**Author Contributions:** All authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been funded by the German Federal Ministry of Education and Research under support code 01LY1218C. The publication of this article was funded by Chemnitz University of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
