Nonlinear dynamics analysis of a low-temperature-differential kinematic Stirling heat engine

The low-temperature-differential (LTD) Stirling heat engine technology constitutes one of the important sustainable energy technologies. The basic question of how the rotational motion of the LTD Stirling heat engine is maintained or lost based on the temperature difference is thus a practically and physically important problem that needs to be clearly understood. Here, we approach this problem by proposing and investigating a minimal nonlinear dynamic model of an LTD kinematic Stirling heat engine. Our model is described as a driven nonlinear pendulum where the motive force is the temperature difference. The rotational state and the stationary state of the engine are described as a stable limit cycle and a stable fixed point of the dynamical equations, respectively. These two states coexist under a sufficient temperature difference, whereas the stable limit cycle does not exist under a temperature difference that is too small. Using a nonlinear bifurcation analysis, we show that the disappearance of the stable limit cycle occurs via a homoclinic bifurcation, with the temperature difference being the bifurcation parameter.

Introduction. -Technology that can use thermal energy as a resource is one of the most important technologies needed to develop a sustainable society. The Stirling heat engine, which was invented in the early 19th century, has received renewed attention over the past several years for its potential application in sustainable energy technology [1,2]. Owing to technological innovations during the past few decades, some types of the Stirling heat engine that can operate under very low temperature differences have been developed [3], and make it possible for the use of thermal energy in our daily life. The demonstration of a low-temperature-differential (LTD) Stirling heat engine operating in the palm of the hand under the difference between the body temperature and room temperature is remarkable [3], and serves as useful material for both teaching and learning thermodynamics.
From a fundamental physics perspective, it is important to understand how the rotational motion of the LTD Stirling heat engine is maintained or lost based on the temperature difference. In contrast to the ideal Stirling thermodynamic cycle with an external operation [1,2], actual Stirling heat engines run autonomously in a selfsustained manner under a sufficient temperature difference. This is accomplished with the aid of a crank-piston mechanism [1,2]. Mathematical modeling of the Stirling the mechanical degrees of freedom of the engine [4]. While a simple theoretical model of an LTD Stirling heat engine that reproduces the rotational motion of the engine has been proposed [5], how this motion is lost with decreasing temperature difference remains to be clarified.
In this study, we provide an answer to this question by proposing and investigating a minimal model of an LTD Stirling heat engine based on nonlinear dynamics. To be more specific, we model a kinematic Stirling heat engine that has been adopted to implement actual LTD Stirling heat engines [2]. Our model is the simplest in the sense that dynamical equations with the fewest degrees of freedom and a stable limit cycle solution interpreted as the rotational motion of the engine are involved. We then analyze a bifurcation of the stable limit cycle solution, with the temperature difference being the bifurcation parameter.
Kinematic Stirling heat engine. -The kinematic Stirling heat engine is a gamma-configuration (cylindersplit) Stirling heat engine [1,2]. In fig. 1, we illustrate a kinematic Stirling heat engine designed for an LTD operation (a prototype of this LTD kinematic engine is the N-92 type designed and built by Senft [3]): The working substance, a gas, is confined into a cylinder with total volume V . The cylinder consists of a large cylinder with volume V d and a small cylinder with volume V p , where V = V d + V p . The small cylinder with a power piston on its top is placed on the large cylinder. The power piston is connected to a crank with a small radius via a connecting rod. The reciprocating motion of the power piston is converted to a rotational motion by the crank. The cylinder is separated into upper and lower regions by a reciprocating displacer piston in the large cylinder. The displacer piston is connected to the other crank with the same radius as the former via a connecting rod. The two cranks, along with a flywheel with a large moment of inertia, are combined to rotate as a single synchronized unit with the same crankshaft centered at the cranks. The phase angle θ d (mod 2π) of the crank connected to the displacer piston is arranged such that it leads the phase angle θ (mod 2π) of the crank connected to the power piston by π 2 as θ d = θ + π 2 [2], where θ d and θ increase clockwise from the origin corresponding to the lowest point of each piston. This unit, which is supported by a pillar, is fixed above the cylinder. There is a small space of negligible volume between the displacer piston and the wall of the large cylinder. The gas is transferred back and forth through the space between the upper and the lower regions of the cylinder because of the reciprocating motion of the displacer piston. The volume V d of the large cylinder that agrees with the swept volume of the displacer piston remains constant, while that of the small cylinder V p changes with θ. At the bottom of the large cylinder, the engine contacts a heat reservoir (e.g., a warm palm or cold ice) at temperature T b . The surrounding air with temperature T air serves as the other heat reservoir and the atmospheric pressure p air acts on the power piston. As shown in fig. 1, an LTD Stirling heat engine has some characteristics in its geometry that include the following [2,3]. The compression ratio of the engine is small; the diameter of the large cylinder and the displacer piston is large; the displacer piston is short, and its stroke is small; the bottom and top surfaces of the large cylinder for thermal conduction is large.
For ∆T ≡ T b − T air > 0, the kinematic Stirling heat engine runs clockwise. Initially, the phase angle of the crank is located at θ = 0 (θ d = π 2 ) with a positive angular velocity. The points θ = 0 and θ = π where the volume of the cylinder reaches its minimum and maximum are called the top dead center (TDC) and the bottom dead center (BDC), respectively [2,6]; a torque force transmitted from the power piston to the crank vanishes at these points. The engine then undergoes the following cyclic process ( fig. 2) [1,2]. (i) → (ii) transfer stroke: The crank rotates clockwise and overcomes θ = 0 (TDC) with the aid of the moment of inertia of the flywheel even against a frictional torque, and the gas expands from the minimum volume V = V min = V d at θ = 0 to push the power piston against the atmospheric pressure. Meanwhile, the displacer piston transfers most of the gas to the lower region of the cylinder. The more gas that gathers in the lower region, the more it is heated up by the hot heat reservoir (warm palm) with T b at the bottom surface of the large cylinder. (ii) → (iii) expansion stroke: The gas temperature becomes close to T b at θ = π 2 (θ d = π) where the gas is completely separated from the cold heat reservoir (surrounding air) at the top surface of the large cylinder. The gas continues to expand and pushes the power piston until the volume reaches the maximum V = V max at θ = π (θ d = 3 2 π). (iii) → (iv) transfer stroke: With the aid of the moment of inertia of the flywheel, the crank keeps rotating, overcoming θ = π (BDC), and the power piston starts to compress the gas against the gas pressure. Meanwhile, the displacer piston transfers the gas to the upper region of the large cylinder. The more gas that gathers in the upper region, the more it is cooled down by the cold heat reservoir with T air at the top surface of the large cylinder. (iv) → (i) compression stroke: The gas temperature becomes close to T air at θ = 3 2 π (θ d = 2π) where the gas is completely separated from the hot heat reservoir at the bottom surface of the large cylinder. The power piston continues to compress the gas until the volume reaches the minimum 2 ) and the engine returns to the initial state. In this cyclic process, the displacer piston serves as an automatic controller that realizes the alternate attachment and detachment of the heat reservoirs with the gas, depending on its location. In contrast, the power piston is a motive part that converts the thermal energy into the rotational motion of the engine via the crank.
This engine runs counterclockwise for ∆T < 0 [3]. The engine starts from V = V min = V d at θ = 0 with a negative angular velocity and follows the following sequence: (i) → (ii) transfer stroke, (ii) → (iii) expansion stroke, (iii) → (iv) transfer stroke, and (iv) → (i) compression stroke, and this sequence is the same as the clockwise case with ∆T > 0 in fig. 2. One difference is that the gas is transferred to the upper region of the cylinder at the first transfer stroke and is heated up by the surrounding air (which acts as the hot heat reservoir) while it is transferred to the lower region of the cylinder at the second transfer stroke and is cooled down by, e.g., cold ice acting as the cold heat reservoir.
Nonlinear dynamic model. -We propose a minimal nonlinear dynamic model that captures the above features of the LTD kinematic Stirling heat engine. We assume that the gas, the working substance, can be approximated as an ideal gas. We also assume that the gas has spatially uniform temperature T and density over the entire cylinder with total volume V .
Our model is a coupled system with mechanical and thermodynamic degrees of freedom [4,5]. The mechanical degree of freedom is the phase angle θ of the crank, which has a radius r. We assume that the equation of motion of the crank is given by where F ≡ σ p nRT V (θ) − p air . Here, I is the moment of inertia of the crank with the flywheel attached, V p (θ) ≡ r(1 − cos θ)σ p and V d = 2rσ d are the volume of the small cylinder and that of the large cylinder, respectively, with σ p and σ d being the bottom areas of these cylinders. The total volume V (θ) of the cylinder is given by where V min = V d = 2rσ d at θ = 0 and V max = V min + 2rσ p at θ = π. F is a net force acting on the power piston created by the pressure difference between the internal gas and the surrounding air. The gas pressure is given by the equation of state for the ideal gas nRT V (θ) with n and R being the number of moles of the gas and the molar gas constant, respectively. Γ is a coefficient of the frictional torque, and the mass of the displacer piston, mass of the connecting rods, and gravity are negligible for simplicity. We can derive eq. (1) from the standard setup of a crank-piston mechanism [6] by imposing conditions r l ≪ 1 and mpr 2 I ≪ 1, where l and m p are the length of the connecting rod and the mass of the power piston with negligible friction, respectively. Under the condition r l ≪ 1, the net force F acting on the power piston is transmitted to the crank via the connecting rod in a nearly parallel fashion [4]. This allows us to interpret the part F sin θ on the r.h.s. of eq. (1) as the tangential component of the force applied to the crank and hence rF sin θ as the rotational torque, which vanishes at θ = 0 (TDC) and θ = π (BDC).
The thermodynamic degree of freedom in our model is the gas temperature T . The time-evolution equation of T is given by the first law of thermodynamics (the law of energy conservation) for the ideal gas [7]: where f ≥ 3 is the total number of spatial and internal degrees of freedom per molecule and dV (θ) We assume that the heat flows from the two heat reservoirs according to the Fourier law and flows only through the top and the bottom surface of the large cylinder into the gas. G j (θ) (j = b, air) in eq. (3) is a statedependent thermal conductance defined by Here G is a usual thermal conductance and is a function expressing a degree of coupling between the gas and the heat reservoir on each side of the large cylinder separated by the displacer piston, where are the separated volumes of the large cylinder by the displacer piston. According to eq. (5), the amount of heat flowing into one side of the large cylinder increases owing to the increase in volume on this side, while the amount p-3 of heat flowing into the other side of the large cylinder reduces with the decrease in volume on that side. The degree-of-coupling function eq. (5) models the role of the displacer piston as shown in the transfer strokes in (i) → (ii) and (iii) → (iv) in fig. 2, where the reciprocating motion of the displacer piston alternately transfers the gas into one side of the large cylinder, making the gas thermally isolated from the heat reservoir on the other side of the large cylinder (χ b ( π 2 ) = 1 and χ air ( π 2 ) = 0 correspond to (ii) in fig. 2 and χ b ( 3 2 π) = 0 and χ air ( 3 2 π) = 1 correspond to (iv) in fig. 2).
To clarify the role of the displacer piston from another perspective, it is convenient to rewrite the heat flow term in eq. (3) as j=b,air where we used eqs. (4)- (6), and is an effective temperature. While this engine runs autonomously, we may interpret eq. (7) as if one heat reservoir with T eff (θ) is tactically attached to the engine by an external agent as in the case of a usual thermodynamic cycle. We can also rewrite eq. (8) as a function of V by using eq. (2) as where the plus (minus) sign applies for 0 ≤ sin θ ≤ 1 (−1 ≤ sin θ < 0) (the solid curve in fig. 3). This is comparable to the T -V diagram of the ideal Stirling thermodynamic cycle [1,2] (the dashed line in fig. 3). An important difference is that the ideal cycle has a square shape with a clear distinction between the temperature-changing process and the volume-changing process whose temperature modulation is a periodic square-shaped function of θ. However, the T eff -V diagram of the present kinematic engine has a circular shape, as in eq. (9), whose temperature modulation is the sinusoidal-shaped function of θ as in eq. (8), which loosely approximates the ideal cycle and is better suited for a practical engine. We note that essentially the same idea of a state-dependent thermal contact of this kind was adopted in earlier studies on different autonomous heat engines [8][9][10].
For convenience, we nondimensionalize our coupled equations (1)  T ≡ T Tair , and ∆T ≡ ∆T Tair is a nondimensionalized temperature difference between the heat reservoirs that will serve as a bifurcation parameter. By defining ω ≡ dθ dt as a nondimensionalized quantity of the angular velocity ω ≡ dθ dt , we can obtain the nondimensionalized dynamical equations corresponding to eqs. (1) and (3) as whereT eff (θ) ≡ 1 + 1+sin θ 2 ∆T andṼ (θ) ≡ 2 +σ(1 − cos θ). For simplicity, we make the assumption of a time-scale separation between θ andω as slow variables andT as a fast variable. Under this assumption,T in eq. (12) starting from an initial value at an initial time quickly relaxes to an adiabatic solutionT (θ,ω) determined by the slow variables, which is obtained by solving dT dt = 0 in eq. (12): By substituting eq. (13) into eq. (11), we can eliminateT from eq. (11) (adiabatic elimination [11]) and obtain the simple dynamical equations closed by the slow variables instead of eqs. (10)-(12) as follows: dω dt =σ (p(θ,ω) −p air ) sin θ −Γω, With eqs. (14) and (15), we can consider the dynamics of our engine on the phase plane (θ,ω) as a dynamical system. This model with eqs. (14) and (15) may be regarded as a certain type of a driven nonlinear pendulum [12] that is motive-powered by a temperature difference.
Numerical calculations. -We perform numerical calculations of eqs. (10)- (12) and eqs. (14) and (15). We choose the parameters in the equations by assuming actual palm-sized LTD Stirling heat engines, some of which have been determined by referring to [3,5]: We set f = 3, r = 0.005 m, I = 0.00025 kg · m 2 , σ d = 0.005 m 2 , σ p = 0.02σ d = 0.0001 m 2 , p air = 1020 hPa, G = 45 J · s −1 · K −1 , T b = 257 ∼ 329 K, T air = 293 K, n = 0.002031 mol, and R = 8.314 J · K −1 · mol −1 . These parameters yield the nondimensionalized parameters in the equations asσ = 0.02,G ≃ 18.94,p air ≃ 0.5154, and ∆T ≃ −0.1229 ∼ 0.1229. Because an estimation of Γ is difficult, here we simply chooseΓ = 0.0025 as an adjusted value that makes the rotational speed of the engine to be approximately dozens to hundreds of revolutions per minute (rpm) on average, as observed in actual LTD Stirling heat engines [3]. We use the 4th Runge-Kutta method with time interval ∆t = 0.0001 in the numerical calculations. Hereafter, we denote by a quantity with · a long-time averaged one.
In fig. 5, we show ω obtained by solving eqs. (14) and (15) as a function of the bifurcation parameter ∆T . The two attractors-the stable limit cycle and the stable fixed point-correspond to the rotational state with ω = 0 and the stationary state with ω = 0, respectively. These states coexist for the sufficiently large |∆T |, while for the too-small |∆T |, the rotational state does not exist. The two bifurcation points where the stable limit cycle disappears are ∆T c ≃ 0.023123 and ∆T ′ c ≃ −0.031305. The diagram is not symmetric as ∆T ′ c = −∆T c because eqs. (14) and (15) are not invariant under the transformation ∆T → −∆T , θ → −θ andω → −ω. For ∆T c < ∆T , the engine rotates clockwise while for ∆T < ∆T ′ c , it rotates counterclockwise, reproducing the actual feature of the engine [3]. While the basin structure yielding the stable limit cycle may be complicated, the initial angular velocities with the sufficiently large magnitude in the expected rotational direction of the engine depending on the sign of ∆T always lead to the rotational states without reaching the stationary states.
Fixed-point analysis. -As a preparation to explain the bifurcations in fig. 5, we investigate fixed points x * ≡ (θ * ,ω * ) of eqs. (14) and (15). There are four total fixed points at most.
As understood from this analysis, at least one stable fixed point always exists, which explains the stationary state with ω = 0 in fig. 5.
Homoclinic bifurcation. -We show that the disappearance of the stable limit cycle at ∆T c and ∆T ′ c in fig. 5 occurs via a homoclinic bifurcation as a global bifurcation [12,13], where the saddle point x * π plays an essential role in this bifurcation.
In fig. 7 (a) and (b), the stable spiral point x * 0 and saddle point x * π are shown on the phase plane for ∆T = 0.025 and ∆T = 0.021 before and after the bifurcation at ∆T c ≃ 0.023123 in fig. 5. x * PE1 and x * PE2 do not exist in these ∆T 's as seen from fig. 6. The stable and the unstable manifolds W s (x * π ) and W u (x * π ) of x * π and the stable limit-cycle orbit with ω > 0 are also shown on the phase plane. Here, W s (x * π ) (W u (x * π )) of x * π is defined as the set of initial values of x(t) satisfying limt →∞ x(t) = x * π p-6 Nonlinear dynamics analysis of a low-temperature-differential kinematic Stirling heat engine (limt →−∞ x(t) = x * π ) [12]. The bifurcation scenario is then described as follows [12]. In fig. 7 (a), we can see that the branch of W u (x * π ) labeled U (the dashed curve) converges to the stable limit-cycle orbit (the bold curve) ast → ∞. As ∆T decreases, these two orbits come close to each other. Finally, at ∆T = ∆T c , the limit-cycle orbit touches the saddle point x * π , forming a homoclinic orbit connecting x * π to itself in an infinite period. At the same time, the limit-cycle orbit disappears. As ∆T further decreases, U becomes a heteroclinic orbit connecting x * π and x * 0 . The same scenario applies for the disappearance of the stable limit-cycle orbit with ω < 0 at ∆T ′ c ≃ −0.031305 in fig. 5. Compare the behavior of U for ∆T = −0.0345 and ∆T = −0.029 before and after the bifurcation in fig. 7 (c) and (d).
This bifurcation is physically understood as follows [12]. As |∆T | decreases, it becomes increasingly difficult for the crank to sustain the rotational motion by overcoming BDC against the frictional torque even with the aid of the moment of inertia. At the critical value, the crank with vanishingly small speed at BDC takes an infinitely long time to pass and then return to it, which corresponds to the formation of the homoclinic orbit. The crank can no longer sustain the rotational motion for the smaller |∆T |.
Discussion. -Our model captures the qualitative features of actual LTD kinematic Stirling heat engines and predicts the homoclinic bifurcation. Accordingly, comparison with experiments must be important to confirm the validity of the model. The bifurcation points ∆T c ≃ 0.023123 and |∆T ′ c | ≃ 0.031305 in fig. 5 correspond to ∆T c ≃ 6.775 K and |∆T ′ c | ≃ 9.172 K in the original parameters. These values less than 10 K may be comparable to measured values at which the actual engines including the N-92 type begin rotating [3]. The model thus reproduces the actual engine behavior quantitatively to some extent, as well as qualitatively. However, it seems that experiments that aim to investigate a type of a bifurcation in LTD kinematic engines have not yet been conducted. Because our bifurcation scenario is a mathematically natural consequence of the physically sound model, it is important to verify this scenario experimentally. Besides, experimental evaluations of the extent to which this minimal model oversimplifies actual engines will also be important, especially for an appropriate extension of the model.
Conclusion. -We proposed a minimal nonlinear dynamic model of the LTD kinematic Stirling heat engine. By replacing the role of the displacer piston with the statedependent thermal conductance and eliminating the thermodynamic degree of freedom adiabatically, we derived our dynamical equations that describe the engine as the driven nonlinear pendulum. Our model reproduced the rotational motion of the engine with the expected direction as the stable limit cycle and the stationary state of the engine as the stable fixed point, respectively. We clarified that the stable limit cycle disappears via the homoclinic bifurcations, with the temperature difference being the bifurcation parameter. Some of the topics that we could not discuss in this study, such as an analysis of the performance of the engine, will be reported elsewhere. Although we have considered the simplest case of the kinematic Stirling heat engine, it is of great interest to investigate a bifurcation scenario of a Ringbom engine [5,14] as the other LTD Stirling heat engine [2] based on nonlinear dynamics. Moreover, we may even invent completely new types of LTD heat engines by using nonlinear dynamics as a powerful design and analytical tool. We expect that the present work will contribute to deepening the understanding of the physics of the Stirling heat engine and facilitate the use of Stirling heat engines for the development of sustainable societies. * * * The author is grateful to H. Kori for helpful discussions. This work was supported by JSPS KAKENHI Grant Number 16K17765.