Global Stability Analysis of a Curzon – Ahlborn Heat Engine under Different Regimes of Performance

We present a global stability analysis of a Curzon–Ahlborn heat engine considering different regimes of performance. The stability theory is used to construct the Lyapunov functions to prove the asymptotic stability behavior around the steady state of internal temperatures. We provide a general analytic procedure for the description of the global stability by considering internal irreversibilities and a linear heat transfer law at the thermal couplings. The conditions of the global stability are explored for three regimes of performance: maximum power (MP ), efficient power (EP ) and the so-called ecological function (EF ). Moreover, the analytical results were corroborated by means of numerical integrations, which fully validate the properties of the global asymptotic stability.


Introduction
Since the pioneering paper published by Curzon-Ahlborn in 1975 [1], most of the studies within the context of finite time thermodynamics (FTT) have been centered on steady-state values of characteristic functions, like power output [2][3][4][5][6][7], efficient power [8,9], ecological function [10] and saving function [11].For instance, the Curzon-Ahlborn (CA) engine [1] (see Figure 1) is characterized by the efficiency η = 1 − (T 2 /T 1 ), with T 1 > T 2 , under a maximum power operating regime [1,7,12,13].However, the energetic description of this kind of model is limited to steady-state situations, with no evident implications for the system's behavior under dynamical changes [7,14,15].In 2001, and attending the need of dynamical models of CA-type engines, Santillán et al. [16] introduced for the first time a local stability analysis for CA-engines, reporting an exponential decaying behavior of perturbations about the steady-state internal temperatures.Later, Guzmán-Vargas et al. [17] investigated the effects of the heat transfer laws and the thermal conductance on the CA-engine's local stability.Páez-Hernández et al. [18] examined the robustness and time delay effects on the performance and stability properties.Barranco-Jiménez et al. [19] studied the local stability of a heat engine, but considering some economic aspects related to the total cost involved in the engine performance.Recently, Reyes-Ramírez et al. [20] investigated the global stability of a CA-engine operating at maximum power output by means of the Lyapunov method.They constructed the Lyapunov function to prove the asymptotic stability behavior about the steady state of the intermediate temperatures in the engine model.
On the other hand, different regimes of performance have been widely studied in the literature in order to characterize the energetic functions of endoreversible engines.Some authors have proposed other maximization criteria, which are quantitatively different compared to the maximum power output; for instance, the minimization of the entropy production, which is represented as the most suitable for the conservation of natural resources.Another regime of performance proposed by Angulo-Brown [10] consists of finding the best trade-off between high power output and low entropy production, this regime is known as the ecological criterion.More recently, Yilmaz [8] proposed a new performance criterion represented by a trade-off between the power output and the efficiency of the system, and it is known as the efficient power regime.
In this work, following the procedure reported in [20], we extend the global stability study considering the other two regimes of performance: maximum efficient power (EP ) [8,9] and maximum ecological function (EF ) [10].The work is organized as follows: In Section 2, we explain the main steady-state characteristics of a heat engine model under different regimes of performance.In Section 3, we describe the global stability analysis method based on Lyapunov's theory to construct the Lyapunov functions for the engine model, both at maximum EP and EF conditions, respectively.In Section 4, we discuss our findings.Finally, in Section 5, we present our conclusions.

Steady-State Characteristics of a Curzon-Ahlborn Engine Model for Different Regimes of Performance
In this section, we present a brief review of the thermodynamical properties of a non-endoreversible Curzon-Ahlborn (CA) heat engine working under different regimes of performance.In Figure 1, a schematic diagram of the irreversible heat engine (CA model) is shown.This engine consists in a Carnot-like thermal engine that works in irreversible cycles and exchanges heat with external thermal reservoirs at temperatures T 1 and T 2 (T 1 > T 2 ).In the steady state, the temperatures of the Carnot-like cycle isothermal branches are x and ȳ; here, overbars are used to indicate the corresponding steady-state value.Under steady-state conditions, heat flows from the hot reservoir to the internal engine and from the engine to the cold reservoir, which are denoted by J1 and J2 , respectively (see Figure 1).Applying the Clausius theorem and using the fact that the inner Carnot-like engine works in irreversible cycles, we get the following inequality, This expression can be transformed into an equality by introducing a parameter R leading to [21,22], The parameter R takes values within the interval 0 < R ≤ 1 (R = 1 corresponds to the endoreversible limit), and it can be seen as a measure of the departure from the endoreversible regime when R < 1 [18,23].If we assume that the heat flows from T 1 to x and from ȳ to T 2 are of the Newton type, then: where κ is the thermal conductance.For simplicity, we have assumed that the heat exchanges take place in conductors with the same thermal conductance κ.The system's steady-state power output and the efficiency can be defined as, and: By using Equations ( 2), ( 3), ( 4) and ( 6), we can write the steady-state temperatures x and ȳ and the power output as [17,20], where τ = T 2 /T 1 .The efficiency that maximizes the power output (Equation ( 9)) is given by [15], It follows that steady-state working temperatures and the power output can be expressed in terms of τ as follows [18], On the other hand, we can construct the expressions of the efficient power [8,9] and the ecological function [10] following a similar procedure to that used to get Equation ( 9).Specifically, these functions are given by PEP = η P and Ē = P − T 2 σ, respectively, where σ is the total entropy production of the cycle.After a little algebra, we arrive at the following expressions, and: where we have applied the second law of thermodynamics to calculate the total entropy production given by σ = J2 T 2 − J1 T 1 (see Figure 1).If we calculate the derivatives of PEP and Ē, with respect to η, and solving for the efficiency, we get, and: Equations ( 16) and ( 17) represent the steady-state efficiencies working under both maximum efficient power (η EP ) ( [8]) and maximum ecological function conditions (η EF )( [10]), respectively.On the other hand, we can use Equations ( 6), (10), ( 16) and ( 17) to find the expression that relates the internal variables x and ȳ to the external temperatures T 1 and T 2 for the M P , EP and EF cases, to get, The role of the internal temperatures is crucial to perform the stability analysis, and for this reason, T 1 , T 2 and P must be written in terms of x and ȳ [16,17].In the case of a maximum efficient power regime, we get, Equivalent expressions were obtained for maximum ecological function regime with the help of a symbolic algebra package, but they are quite lengthy to be included.

Dynamical CA-Engine Model
The dynamical model of the CA-engine is based on the rate of change of the internal temperatures, x and y (see Figure 1), that is these temperatures are allowed to change with time, according to the amount of heat entering or exiting [16].The model assumes that now x and y are considered macroscopic bodies with finite heat capacity C. Thus, the following coupled differential system represents the change of the internal temperatures [16], and, where J 1 and J 2 represent the heat flows from the hot reservoir to the internal engine and from the engine to the cold reservoir, respectively.Moreover, assuming that the engine is internally reversible (endoreversible hypothesis), the fluxes J 1 and J 2 , are given by, and, where P represents the power at the regime of performance, with i = M P , EP and EF .We assume that P depends, out of the steady state, on x, y temperatures in a similar way as it depends on x, ȳ.

Lyapunov Functions of the CA-Engine Model under Both Maximum Efficient Power and Maximum Ecological Function
In general, there is not a defined method to find a Lyapunov function, however, a candidate Lyapunov function must satisfy, (i) V (x, y) must be positive definite in a region around the steady state, (ii) V (x, y) must be negative definite, that is, V (x, y) ≤ 0.
Here, we apply Krasovskii's method [24] to construct the mentioned Lyapunov functions.This method consists of finding a definite-negative matrix of the form J(X) = A(X) + A + (X), with X = (x, y) and:  and (x, y) is the fixed point of the system.The candidate Lyapunov function is given by [24], (26)

Maximum Efficient Power Regime
By substituting Equations ( 21), ( 24) and ( 25) into Equations ( 22) and ( 23), we get the following dynamic equations for the temperatures x and y under EP conditions, and: These equation systems are of the form ẋ = f (x, y), ẏ = g(x, y), and have a fixed point at (x,ȳ), that is, f (x, ȳ) = 0 and g(x, ȳ) = 0. Following a similar procedure of [20], we construct the matrix J(X), Due to the fact that the matrix J EP is definite negative, as we can easily verify, the candidate Lyapunov function for the CA-engine model at maximum efficient power is of the form, From the previous equation, as in the case of maximum power conditions [20], we can verify that V EP (x, ȳ) = 0, that is the steady-state regime is the extremum and global minimum of V EP (x, y).Besides, the previous equation satisfies,

Maximum Ecological Function Regime
Following a similar procedure as in the previous section, the dynamic equations for x and y, but now at maximum EF conditions, are given by, and: The Lyapunov function under maximum ecological function conditions is given as, where Λ(x, y, R) = √ R(Rx 2 +8y 2 )−R(x+2y) Ry . For this case, as occurred with Equation (30), it is easy to verify that VEF (x, y) = ∂V EF ∂x ẋ + ∂V EF ∂y ẏ = ∇V EF , Ẋ ≤ 0.

Discussion
We have showed that the functions V EP (x, y) and V EF (x, y), given by Equations ( 29) and (33), which correspond to both maximum EP and EF conditions, respectively, are the Lyapunov functions for their corresponding dynamical system, that is, in both cases, the system is asymptotically stable [24].To further explore the form of the functions V EP (x, y) and V EF (x, y), a set of curves are constructed to illustrate their behavior for the feasible region.Figures 2 and 3 show the level curves defined by V EP (x, y) = α EP and V EF (x, y) = α EF , respectively, for different values of the constant α EP and α EF .As expected, we observe in both cases that, as the constant value decreases, the surface also decreases towards the steady-state values (x, ȳ), for the endoreversible and non-endoreversible cases, respectively.Besides, in Figures 4 and 5, we show qualitative surface plots of V EP (x, y) and V EF (x, y).As we can see, for both cases, the surfaces resemble a basin, which indicates that the global stability of the system is preserved for different values of the irreversibility parameter R under both maximum EP and maximum EF conditions.Additionally, we performed numerical integrations of the system to further verify the global asymptotically properties for the regimes.In Figure 6, we present the results of numerical integrations of the system given by Equations ( 27) and (28) for the case of efficient power and Equations ( 31) and (32) for the case of ecological function.For illustrative purposes, we have considered four representative cases for the initial conditions (see the caption of Figure 6 for details).As we can see, in all cases, the internal temperatures evolve towards the steady state.For convenience, in this figure, we also include the case of M P previously reported in [20].In order to compare the decaying rates of the internal temperatures, we also calculate the rate of change, dx dt and dy dt , for the evolution curves shown in the main frame of Figure 6.The results are depicted in the insets of Figure 6.We observe that for the Cases (1) and ( 2), the rate of change of x and y corresponding to the M P -regime are bigger than the corresponding to EF and EP , indicating that a fast approximation takes place under M P conditions, whereas for Cases (3) and ( 4), the opposite situation is observed.We recall that these results are valid for the particular values considered to perform the numerical integrations.We have considered four representative cases, where each has different initial conditions (x i 0 , y i 0 ) (which will be useful for the description of Figure 7), with i = 1, 2, 3, 4. In all of the cases, both temperatures exhibit a decaying behavior towards their corresponding steady-state values, and we observe that x(t) approximates faster than y(t) for the four cases.To obtain the evolution as a function of time, the equations were numerically integrated using the Runge-Kutta triple (Bogacki et al. [25,26]).(Insets) First derivative with respect to time of the curves in the main frame.The rate of change reveals that fast approximations towards the steady state are observed for MP-conditions (Cases (1) and ( 2)), while the opposite situation is valid for Cases (3) and (4).

R=1
case (1) case (2) case ( 3) case (4) x 0 (2) x 0 (3) x 0 (4) y 0 (4) y 0 (1) y 0 (2) Moreover, we construct phase spaces for the cases R = 1, R = 0.8 and R = 0.6 to observe the trajectories' approximation towards the steady state, that is to observe the way internal temperatures approximate their corresponding steady sate for each regime (Figure 7a-c).In all cases, the temperatures display a fast rate along the horizontal axes in comparison with the vertical approximation.Furthermore, from Figure 7a-c, it is noteworthy that the location of the steady state moves to the right and slightly downward as R decreases and also when the operation regime changes in the order M P , EP and EF .We notice that these approximations observed in the phase portrait are related to the existence of fast and slow eigendirections when a local stability analysis is performed around the steady-state values [17,20].
Figure 7. Qualitative phase portrait of x(t) vs. y(t) for the endoreversible case R = 1 (a) and the non-endoreversible cases R = 0.8 (b) and R = 0.6 (c) considering the three regimes performance: M P , maximum EP and maximum EF conditions.In (a), the trajectories represent the four cases described in Figure 6.We use the same initial points to show the evolution towards the steady-state values.Notice that the steady-state coordinates move to the right and downward of their location for each regime and for the different values of R. x

Concluding Remarks
The global stability of a CA-engine was analyzed by considering three regimes of performance: maximum power, efficient power and ecological function.The Lyapunov analysis presented here represents an attempt to characterize the global properties of the CA-engine model.Our results show that the stability is preserved for a range of parameter values, and the functions V EP (x, y) and V E (x, y) do not change its stability properties, even in the case when irreversibilities are present, that is for R < 1.Moreover, the findings of the global asymptotically-stable behavior of the equilibrium state were corroborated by means of numerical integrations performed over the dynamical equations of the system.The internal temperatures approximate the steady-state values, confirming that the system is stable for different initial configurations.Our results confirm the global stability properties of the CA-engine for three regimes of performance.
matrix (A + represents the transpose matrix) and where F (x, y) = (f, g),

Figure 4 .
Figure 4. Qualitative surface plots of the Lyapunov function V EP (x, y), for different values of the irreversibility parameter R. We observe that in all cases, the surfaces are a basin shape-type, even when R takes small values (R < 1).

Figure 5 .
Figure 5. Qualitative surface plots of the Lyapunov function V E (x, y), for different values of the irreversibility parameter R. As in the case of V EP (x, y) (Figure 4), surfaces also are the basin shape-type.

Figure 6 .
Figure 6.(Main frame) Evolution of working temperatures vs. time t, for different initial conditions for the endoreversible case R = 1 (the cases R = 0.8 and R = 0.6 are not shown).We have considered four representative cases, where each has different initial conditions (x i 0 , y i 0 ) (which will be useful for the description of Figure7), with i = 1, 2, 3, 4. In all of the cases, both temperatures exhibit a decaying behavior towards their corresponding steady-state values, and we observe that x(t) approximates faster than y(t) for the four cases.To obtain the evolution as a function of time, the equations were numerically integrated using the Runge-Kutta triple (Bogacki et al.[25,26]).(Insets) First derivative with respect to time of the curves in the main frame.The rate of change reveals that fast approximations towards the steady state are observed for MP-conditions (Cases (1) and (2)), while the opposite situation is valid for Cases (3) and (4).