Synchronization approach to achieving maximum power and thermal efficiency for weakly-coupled low-temperature-differential Stirling engines

Low-temperature-differential (LTD) Stirling engines are heat engines that can operate autonomously with a slight temperature difference between low-temperature heat reservoirs and are thus expected to contribute to a sustainable society. A minimal dynamical-system model with only two variables has been proposed to explain the principle of autonomous rotational motion caused by temperature differences, and the maximum efficiency of the engine was formulated [Y. Izumida, Europhys. Lett. 121, 50004 (2018); Phys. Rev. E 102, 012142 (2020)]. This paper aims to clarify the coupling effects on the dynamics, power, and thermal efficiency of a pair of weakly coupled LTD Stirling engines and formulate the maximum thermal efficiency of the coupled system in the quasilinear response regime. We show that the dependence relation between the effective frequency difference and the coupling strength is characterized by a hysteresis, which comes from different kinds of bifurcations in the process of increasing and decreasing the value of the coupling strength. Then, by generalizing thermodynamic fluxes and forces and their quasilinear relations for engines under weak coupling, we show that the coupling improves the power exerted against the load torques and the thermal efficiency. We further show that their maximum values are achieved when the engines are synchronized. Since the thermal efficiency depends on the frequency difference, the dependence of thermal efficiency on the coupling strength is also characterized by a hysteresis. Finally, the load torque that achieves the maximum thermal efficiency of the coupled system is formulated.


I. INTRODUCTION
A heat engine is a system that uses thermal energy from a high-temperature heat reservoir to extract positive work.According to the second law of thermodynamics, a low-temperature heat source is required to discard part of the thermal energy to extract positive work from a heat engine.Low-temperature differential (LTD) Stirling engines, which can rotate autonomously with only a slight temperature difference between low-temperature heat reservoirs, are attracting significant attention as an elemental technology to realize a sustainable society [1][2][3].From this perspective, it is vital to understand the dynamical characteristics of LTD Stirling engines through appropriate mathematical modeling [4,5].A nonlinear dynamics model has been proposed to explain the loss of rotational motion of LTD Stirling engines, which was found to be caused by a homoclinic bifurcation [6].
Another important issue for the LTD Stirling engines is thermal efficiency.In [7], one of the authors demonstrated that the engine's rotational state is in a quasilinear response regime where the thermodynamic fluxes show a linear dependence on the thermodynamic forces and formulated the maximum efficiency of the engine based on the fact that the response coefficients of the quasilinear relations are symmetric, which is similar to Onsager symmetry in linear irreversible thermodynamics.However, the power extracted from a single LTD Stirling engine is quite limited, thus it is desirable to operate a population of Stirling engines to extract adequate work for practical purposes.Methods that achieve maximum efficiency by properly controlling a population of Stirling engines then turn out to be important.
Synchronization is a self-organized phenomenon in which oscillators align their rhythms through interaction and is widely observed in natural and artificial systems [8,9].A natural question that would be raised is whether synchronization through coupling between the LTD Stirling engines can improve the total power and thermal efficiency.If it does, then higher power and thermal efficiency can be achieved simply by allowing the engines to interact with each other.Although experimental studies on the synchronization of LTD Stirling engines have been conducted [10,11], the effects of synchronization on power and thermal efficiency have not yet been clarified theoretically.
This paper aims to clarify the coupling effects on the dynamics, power, and thermal efficiency of a pair of weakly coupled LTD Stirling engines and to formulate the maximum thermal efficiency of the coupled system in the quasilinear response regime.We will provide a model of a pair of weakly coupled LTD Stirling engines and investigate the coupling effects on the dynamics through numerical experiments.After that, we will provide a theoretical analysis of the effects of the weak coupling on power and thermal efficiency.By generalizing thermodynamic fluxes and forces and their quasilinear relations for engines under weak coupling, we show that the coupling improves the power exerted against the load torques and the thermal efficiency.We further show that their maximum values are achieved when the engines are synchronized.Finally, we formulate the load torques that achieve the maximum thermal efficiency of the coupled system.

II. MODEL
We consider a pair of weakly coupled LTD Stirling engines with the same parameters except for the load torques T (1) load and T (2) load acting on the cranks (Fig. 1).Heat reservoirs at temperatures Tb and Tt ( Tb > Tt ) are attached to the bottom and top surfaces of the large cylinders of the engines respectively, and we define the temperature difference ∆ T ≡ Tb − Tt for later use.The temperature difference ∆ T and load torque T (i) load (i = 1, 2) are assumed to be sufficiently small.A nondimensionalized minimal model of a single LTD Stirling engine has been proposed in [6] with the following form: where θ is the phase angle of the crank connected to the power piston; σ is a positive constant determined by the surface areas of the large and small cylinders; represent the volume and temperature of the gas confined to the cylinders, respectively; T eff (θ) = 1 + sin θ 2 ∆T is the effective temperature of the heat reservoirs that periodically changes depending on the phase angle; G is the thermal conductance associated with the heat transfer between the gas and the surface of the large cylinder; P air is the atmospheric pressure acting on the power piston, and Γ is the friction coefficient associated with the power piston.All the variables and parameters that are off-tilde represent dimensionless quantities.The minimal model was obtained by assuming that the heat fluxes from the bottom and top surfaces of the large cylinder obey the Fourier law , where G m (θ) with m = b (or t) represents the effective thermal conductance between the gas and the bottom (or top) heat reservoir.It was also assumed that G m (θ) ≡ Gχ m (θ), where χ m (θ) (0 ≤ χ m (θ) ≤ 1) is a function that controls the coupling between the gas and the bottom or top heat reservoir, given as χ b (θ) = 1 2 (1 + sin θ) and χ t (θ) = 1 2 (1 − sin θ) [7].The dynamical equations describe the engines as coupled nonlinear pendulums, where the first term on the RHS of Eq. (1b) represents the driven force due to the temperature difference.Since it has been experimentally demonstrated that the minimal model (1a)-(1b) explains the essential characteristics of a real LTD Stirling engine [12], we generalize the above minimal model by adding a coupling term to describe the dynamics of a pair of weakly-coupled LTD Stirling engines i and j (i, j ∈ {1, 2}, i j): The last term in Eq. (2b) represents the coupling with K > 0 being the coupling strength.Note that the coupling should be anti-symmetric according to the action-reaction law and is chosen to be a sine function for simplicity.

III. COUPLING EFFECTS ON THE DYNAMICS
To evaluate the degree of synchronization caused by the coupling, we introduce the effective frequency as where ... ≡ lim τ→0 1 τ τ 0 ...dt denotes a long-time average and is reduced to the average over one period for engines in periodic motion.For K = 0, the engines are adjusted to be in the quasilinear response regime [7] so that they rotate autonomously in a self-sustained manner.The phase space of engine i is a set of ordered pairs {(θ i , ω i ) : θ i ∈ [−π, π), ω i ∈ R} which is a one-dimensional cylinder T × R, and the rotational motion is described by a limit cycle that circles the surface of the cylinder.The limit cycle for K = 0 is referred to as the unperturbed limit cycle.Since each engine behaves as a limit-cycle oscillator, we can define the natural frequencies of the two engines and denote them by ω (i)  n where i = 1, 2, i.e., ω (i) n = ω i for K = 0. We assume that the natural frequency difference ∆ω n ≡ |ω (1)  n − ω (2)  n | is sufficiently small compared to ω (i) n , but sufficiently larger than the ω i -directional amplitude of the unperturbed limit cycle in the phase space, i.e., ω var i ≪ ∆ω n ≪ ω (i)  n , where ω var i is given by and Ω UNP i denote the set of ω i -components of the points on the unperturbed limit cycle of engine i.For K > 0, the two engines are coupled with each other and synchronization occurs for a sufficiently large coupling strength.
Figure 2 shows the dependence relation between the effective frequency difference ω d ≡ ω 1 − ω 2 and the coupling strength K, where the forward (backward) process corresponds to the situation in which the value of K is increased (decreased).Typical trajectories in the (θ i , ω i ) plane are also shown in the same figure.
The bifurcation diagram in Fig. 2 reminds us of the dynamics of a driven pendulum [13] or a two-node power grid model consisting of one generator and one consumer [14,15], where a homoclinic and saddle-node bifurcation for fixed points occur in the forward and backward processes respectively.Our numerical analysis indicates that similar bifurcations occur in weakly coupled LTD Stirling engines.Particularly, when considering the differential system of (2a)-(2b), a homoclinic bifurcation due to the annihilation of a quasi-periodic trajectory is thought to occur in the forward process as a result of the collision of this quasi-periodic trajectory with a saddle limit cycle corresponding to an unstable synchronous state, and a saddle-node bifurcation is thought to occur in the backward process due to the collision of a stable limit cycle corresponding to a stable synchronous state and an unstable limit cycle corresponding to an unstable synchronous state (See Appendix A for details).

IV. COUPLING EFFECTS ON THE POWER AND THERMAL EFFICIENCY
The power and thermal efficiency of a single LTD Stirling engine in a quasilinear response regime have been derived in [7].Before discussing the coupling effects on the power and thermal efficiency of the total system, we generalize the thermodynamic fluxes and forces as well as their quasilinear relations for engines under weak coupling.
The instantaneous power P (i) produced by engine i is given by where d dt 1 2 ω 2 i is the change rate of rotational energy, P air dV dt is the power that is carried out against the atmospheric pressure, P (i) load ≡ T (i) load ω i is the power that is carried out against the load torque, P (i) fric ≡ Γω 2 i is the power that is carried out against the friction torque, and P (i) K ≡ K sin(θ i − θ j )ω i is the power due to the weak coupling.The power that is carried out against the load torque P (i)  load is referred to as the brake power [16] made by engine i. = 0, the timeaveraged power of engine i is obtained as P (i) = P (i) load + P (i) fric + P (i) K .It should be noted that for the system in quasiperiodic motion, the trajectory is not closed, so the long-time average can not be reduced to the average over an oscillation period.Given that the coupling is sufficiently weak and each engine is in the quasilinear response regime when there is no coupling, ω i can be approximated to be the effective frequency ω i when considering P (i) K , i.e., P (i) K can be approximated as P (i) K ≈ K sin(θ i − θ j ) ω i .Since the time-averaged changes in the entropy and the internal energy of the gas confined to the cylinder are zero, the time-averaged entropy production rate of the total thermodynamic system dσ dt is the sum of the load = 8.5324×10 −7 , and T (2)  load = 1.2799×10 −5 .The values of the parameters other than the load torques are set as the same as Ref. [7], corresponding to the situation where a standard LTD Stirling engine is placed at a temperature difference of a few degrees.The load torques T (1)  load and T (2)  load are chosen so that ω var i ≪ ∆ω n ≪ ω (i) n .Similar graphs can be obtained with a different set of parameter values that satisfy the above conditions.
time-averaged entropy change rates of the two heat baths, which is calculated as = ω m −T (1)  load − T (2)  load where we have used the energy conservation law load + P (i) fric + P (i) K in Eq. ( 6) and approximated T b and T t as their mean value in Eq. ( 7), which equals 1 for the nondimensionalized case.Here, ω m ≡ 1 2 ( ω 1 + ω 2 ) is the mean effective frequency, ω d = ω 1 − ω 2 is the effective frequency difference, and J (1)  Q b + J (2) Q b is the total heat flux from the high-temperature heat reservoir.
Equation (8) suggests that −T (1)  load − T (2) load , −K sin(θ 1 − θ 2 ) − 1 2 T (1) load − T (2)  load , and ∆T can be considered as thermodynamic forces with conjugate fluxes ω m , ω d , and J (1)  Q b + J (2) Q b under appropriate conditions, for which the quasilinear relations are obtained as follows (See Appendix B for details): Here, L 11 , L 12 , L 21 , and L 22 correspond to the quasilinear response coefficients of a single engine in the non-coupling case [7]: where ... θ ≡ 1 2π 2π 0 ...dθ denotes a phase average.We now consider the coupling effects on the averaged brake power P load ≡ T (1)  load ω 1 + T (2)  load ω 2 and thermal efficiency η ≡ Q b by using the generalized quasilinear relations (9) between thermodynamic fluxes and forces.To that end, we rewrite the averaged brake power P load in the following form: P load = T (1)  load ω m + Here, P m ≡ ω m T (1) load + T (2)  load denotes the power owing to the motion of the mean angle and P rel ≡ 1 2 ω d T (1)  load − T (2)   load denotes the power owing to the relative motion.Since ω m is independent of the coupling from Eq. ( 9), we need only consider the coupling effects on P rel .Without loss of generality, we assume T (1) load < T (2)  load , in which case the value of ω d decreases due to the effect of the coupling strength K in both forward and backward processes, as was shown in Fig. 2.This leads to the fact that P rel is an increasing function of K, which means that the coupling improves the averaged power.To see the coupling effect on the thermal efficiency, we notice that the total heat flux from the high-temperature heat reservoir J (1) Q b is independent of the coupling from Eq. ( 9).This suggests that the coupling improves both the averaged brake power and the thermal efficiency given that different load torques act on the cranks, and their maximum values are achieved when the engines are synchronized.
To give a physical interpretation of the fact that a weak coupling improves P rel , let us concentrate on ω d in Eq. ( 9).We find that ω d in the non-coupling case is reduced to −L 11 T (1)  load − T (2)  load , which means that P rel is generated by the synergy of the load torque difference and the relative motion due to the load torque difference when there is no coupling.In this case, P rel takes a negative value as long as T (1)   load T (2)  load , and is a decreasing function of |T (1)  load − T (2)  load |.We thus conclude that P rel reduces the averaged brake power given a fixed sum of load torques.When there is a coupling added, P rel is obtained by Here, ∆P rel ≡ −L 11 K sin(θ i − θ j ) T (1)  load − T (2)  load represents the change of P rel due to the coupling, which takes a positive value as long as T (1)   load T (2)  load .This suggests that the coupling improves the averaged brake power.From Eq. ( 16), we find that the increase in averaged brake power is due to the suppression effect of coupling on relative motion caused by the load torque difference.The averaged brake power takes the maximum value when K sin(θ i − θ j ) reaches − 1 2 T (1) load − T (2)  load , in which case ω d = 0, meaning that the engines are synchronized.
Figure 3 shows the dependence relation between the thermal efficiency and the coupling strength in the forward and backward processes.The blue line is obtained by numerical experiment, while the orange line is obtained by approximate calculation using the quasilinear relations between thermodynamics fluxes and forces.Since it is difficult to calculate K sin(θ 1 − θ 2 ) analytically, we used numerical values of it in the approximate calculation.We can see some gap between experimental and theoretical values, which is caused by neglecting higher-order terms and by the averaging approximation made in the derivation of Eq. ( 9).We also find that the dependence of the thermal efficiency on the coupling strength is characterized by a hysteresis as in the case of the frequency difference in Fig. 2 (a).This is because the thermal efficiency depends on the effective frequencies of the two engines.Such a hysteresis structure facilitates the robustness of maintaining maximum thermal efficiency.K is increased from 0 to 6.0 × 10 −5 in increments of 4.0 × 10 −7 in the forward process and decreased in the same way in the backward process.In the current case, we can confirm that the coupling has increased the thermal efficiency of the total system by about 7%.
We have confirmed that coupling can improve the averaged brake power and thermal efficiency.Since the total load torque determines the power and thermal efficiency given a fixed coupling strength, it is important to investigate the total load torque that achieves the maximum values of them for synchronized engines.In this case, ω m and ω d are reduced to the synchronized frequency ω s , and 0, respectively, indicating that ω s and J (1)  Q b + J (2) Q b are the only thermodynamic fluxes for the coupled system with conjugate forces −T (1)  load − T (2)  load and ∆T .The thermal efficiency is given by which is completely determined by the thermodynamic fluxes and forces.Therefore, the formulation of the maximum thermal efficiency of a single engine given in [7] is directly applicable to the present case.The maximum thermal efficiency η max and the total load torque T (1)  load + T (2)  load that achieves this maximum thermal efficiency are given by where is the coupling-strength parameter and η C ≡ 1 − T t T b is the Carnot efficiency, i.e., the maximum thermal efficiency that a heat engine may have operating between two heat reservoirs.We find that the coupling-strength parameter, as well as the maximum thermal efficiency, is of the same form as that of a single engine, while the total load torque that achieves the maximum thermal efficiency is twice as large as that of a single engine given a fixed temperature difference ∆T .The load torques achieving the maximum power and the corresponding thermal efficiency [17] can be discussed in the same way [7].

V. DISCUSSION AND CONCLUSIONS
In this paper, we have considered a minimal dynamical-system model of weakly coupled LTD Stirling engines and analyzed the coupling effects on the dynamics, power, and thermal efficiency.We clarified the mechanism of different kinds of bifurcation in the forward and backward processes and generalized the thermodynamic fluxes and forces and their quasilinear relations when the weak coupling is incorporated.Based on the linear relations, we concluded that the coupling improves the power exerted against the load torque as well as the thermal efficiency and that their maximum values are achieved when two engines are synchronized.We formulated the maximum thermal efficiency given that the coupled engines are synchronized and found that the expression of the maximum thermal efficiency is given in the same form as that of a single engine.Although the thermal efficiency of LTD Stirling engines is low [18], their great value lies in their ability to generate power with only a small temperature difference.In other words, unlike the large engines that are run in factories, they do not need fuel to generate power, and only a ubiquitous temperature difference (e.g., between air and ground) is needed for the engine to generate power.To achieve sufficient power for practical use, it is desirable to operate a large number of LTD Stirling engines, and synchronizing the engines may be an effective way to further improve power and thermal efficiency.This study discusses the effects of synchronizing two engines as the simplest case, but will be extended to the case of multiple (3 or more) engines in the future.
To gain more insight into the bifurcations that occur when changing the coupling strength, we plot the trajectories in the subspace {(θ d , ω d )} in the forward and backward processes before and after the bifurcations occur, where θ d ≡ θ 1 − θ 2 and ω d ≡ ω 1 − ω 2 .In the forward process, when K is raised to a value slightly less than the bifurcation point K fd , the quasi-periodic trajectory evolves much more slowly near θ d = 3 than elsewhere (See Fig. 4. (a)); After the bifurcation, the trajectory converges to a stable limit cycle (See Fig. 4. (b) and (c)).In the backward process, the stable limit cycle does not disappear until K reaches another bifurcation point K bd ; After the bifurcation, the stable limit cycle collapse and the trajectory converges to a quasi-periodic attractor circling the phase cylinder.These results suggest that a homoclinic bifurcation and a saddle-node bifurcation occurs in the forward and backward process respectively: in the forward process, the quasi-periodic trajectory evolves in the neighborhood of the stable and unstable manifolds of a saddle limit cycle corresponding to an unstable synchronous state when K is slightly smaller than K fd , and converges to a stable limit cycle corresponding to a stable synchronous state after the bifurcation occurs; in the backward process, a saddle-node bifurcation occurs due to the collision of the stable and unstable synchronous states.

APPENDIX B: DERIVATION OF EQ. (9)
We derive the quasilinear relations between thermodynamic fluxes and forces given by Eq. ( 9).To that end, we average both sides of Eq. (2b): where we have used the fact that The trajectory evolves in the neighborhood of a homoclinic orbit starting and ending at a saddle limit cycle corresponding to an unstable synchronous state before converging to a stable limit cycle corresponding to a stable synchronous state.(c) Enlarged view of the stable synchronous state.
since the trajectory remains in the neighborhood of the unperturbed limit cycle.By expanding T (θ i , ω i ) w.r.t.ω i as the first term on the RHS of Eq. ( 21) can be obtained as σ T (θ i , ω i ) V(θ i ) − P air sin θ i = σ T eff (θ i ) V(θ i ) − σ sin θ i GV 2 (θ i ) ω i − P air sin θ i + O(∆T ω i , ω 2 i ) .

FIG. 1 .
FIG. 1.(a) Front view of an LTD Stirling engine.(b) Side view of a pair of weakly coupled LTD Stirling engines with different load torques acting on the cranks.The gases confined to the cylinders are in contact with the bottom and top heat reservoirs.

FIG. 2 . 1 V( π 4 )
FIG. 2. (a)Dependence relation between the effective frequency difference ω d and coupling strength K.In the forward (backward) processes, K is increased (decreased) between 0 and 6.0×10 −5 with a step of 4.0×10 −7 .We can confirm that the transition point K fd in the forward process is different from the transition point K bd in the backward process.Typical trajectories in the forward process are shown in the (θ i , ω i ) plane: (b) K = 0.The trajectory of each engine is a limit cycle that circles the phase cylinder.(c) K = 1.60 × 10 −5 .The trajectory of each engine evolves quasi-periodically and vibrates repeatedly in the vertical direction around a certain value.(d) K = 4.1522 × 10 −5 .The trajectory of each engine is a limit cycle whose periods are the same, i.e., ω 1 = ω 2 , indicating that the two engines are synchronized.Other parameters are set as follows in all subsequent numerical experiments: σ = 0.02, p air = 1V( π 4 )≈ 0.49854, G = 1.5, Γ = 0.001, ∆T = 1/29.3,T (1) load = 8.5324×10 −7 , and T(2)  load = 1.2799×10 −5 .The values of the parameters other than the load torques are set as the same as Ref.[7], corresponding to the situation where a standard LTD Stirling engine is placed at a temperature difference of a few degrees.The load torques T(1)  load and T(2)  load are chosen so that ω var i ≪ ∆ω n ≪ ω (i) n .Similar graphs can be obtained with a different set of parameter values that satisfy the above conditions.

FIG. 3 .
FIG.3.Dependence relation between the thermal efficiency and the coupling strength for (a) forward process and (b) backward process.K is increased from 0 to 6.0 × 10 −5 in increments of 4.0 × 10 −7 in the forward process and decreased in the same way in the backward process.In the current case, we can confirm that the coupling has increased the thermal efficiency of the total system by about 7%.

FIG. 4 .
FIG.4.Trajectories in the subspace {(θ d , ω d )} in the forward process for different values of K. (a) K is set slightly smaller than K fd and the quasi-linear trajectory evolves much more slowly near θ d = 3 than elsewhere.(b) K is set slightly larger than K fd .The trajectory evolves in the neighborhood of a homoclinic orbit starting and ending at a saddle limit cycle corresponding to an unstable synchronous state before converging to a stable limit cycle corresponding to a stable synchronous state.(c) Enlarged view of the stable synchronous state.