*Article* **Optimization, Stability, and Entropy in Endoreversible Heat Engines**

**Julian Gonzalez-Ayala 1,2,\*, José Miguel Mateos Roco 1,2, Alejandro Medina 1,2 and Antonio Calvo Hernández 1,2**


Received: 28 September 2020; Accepted: 18 November 2020; Published: 20 November 2020

**Abstract:** The stability of endoreversible heat engines has been extensively studied in the literature. In this paper, an alternative dynamic equations system was obtained by using restitution forces that bring the system back to the stationary state. The departing point is the assumption that the system has a stationary fixed point, along with a Taylor expansion in the first order of the input/output heat fluxes, without further specifications regarding the properties of the working fluid or the heat device specifications. Specific cases of the Newton and the phenomenological heat transfer laws in a Carnot-like heat engine model were analyzed. It was shown that the evolution of the trajectories toward the stationary state have relevant consequences on the performance of the system. A major role was played by the symmetries/asymmetries of the conductance ratio *σhc* of the heat transfer law associated with the input/output heat exchanges. Accordingly, three main behaviors were observed: (1) For small *σhc* values, the thermodynamic trajectories evolved near the endoreversible limit, improving the efficiency and power output values with a decrease in entropy generation; (2) for large *σhc* values, the thermodynamic trajectories evolved either near the Pareto front or near the endoreversible limit, and in both cases, they improved the efficiency and power values with a decrease in entropy generation; (3) for the symmetric case (*σhc* = 1), the trajectories evolved either with increasing entropy generation tending toward the Pareto front or with a decrease in entropy generation tending toward the endoreversible limit. Moreover, it was shown that the total entropy generation can define a time scale for both the operation cycle time and the relaxation characteristic time.

**Keywords:** multiobjective optimization; Pareto front; stability; maximum power regime; entropy behavior

## **1. Introduction**

The optimization of energy converters has never been as relevant as it is now. Energy production requirements, efficient use of heat sources, and heat waste reduction are continually pushing technological edges. Linked to this degree of specialization are the control and stability of operation regimes yielding a desirable stable production, despite the possible variations of external or internal conditions. In many cases, this requires the fine-tuning of control parameters with intrinsic energetic costs. In this regard, there are some hints as to the possibility of seizing the stability of heat engine operation regimes to enhance their performance by relaxing the control over the operation parameters [1–3]. Studies regarding the weakly dissipative limit [4] thus far have pointed out the special role played by the endoreversible model [5,6]. These studies consider that after the system experiences a perturbation in its operation variables, the trajectories that lead the system back to its steady state tend to evolve toward the endoreversible limit. This limit has been proven to be associated with thermodynamic states, with the best compromise between maximum efficiency, maximum power output, and minimum entropy production as given by the Pareto front of the system. Although the equivalence of the low dissipation model (based on entropy generation) and the endoreversible model (based on specific heat transfer laws) has been established for several heat transfer laws [7–9], it is not obvious that the endoreversible behavior appears interconnected to the stability of the low dissipation model. A natural concern is to look for this very same behavior in the context of endoreversible heat engines.

On the other side, the stability of heat engines is not a novel topic. Especially in the realm of finite-time thermodynamics, endoreversible and irreversible models have, in their assets, a good number of works in this regard. From the pioneering work of Santillan et al. [10], a number of studies have analyzed the local and global stability of a variety of operation regimes [11–18], including economic factors [19], and have extended the analysis to heat pumps, refrigerators, and generalized heat engines [17,20–27].

The validation and applicability of the endoreversible hypothesis has been widely analyzed and discussed in the specialized literature [28–31], and although it constitutes an idealization of an irreversible device, recent studies on molecular dynamics simulation [32] have validated its predictions for a finite-time Carnot cycle of a weakly interacting gas (considered as a nearly ideal gas) and considering the local equilibrium for a Maxwell–Boltzmann distribution with a spatially uniform temperature and a spatially varying local center-of-mass velocity. In particular, these results also point to the validity of the paradigmatic Curzon–Ahlborn efficiency at maximum power [33,34]. These molecular dynamics simulations of a two-dimensional Carnot engine allowed to investigate not only the optimization of power output, but also some other figures of merit involving entropy production as the ecological figure of merit [32]. This reinforces the validity of the Carnot-like endoreversible model, where the quasistatic conditions linked to endoreversibility rely on the thermalization due to the internal dynamic speed. Lastly, a recent work reported some new conceptual insights (simulation and reconstruction) on the endoreversibility hypothesis [35] by considering that subsystems are out of equilibrium, i.e., including internal irreversibilities: First, by means of the contact temperature of the heat flows and the non-equilibrium molar entropy for material flows. This new feature is beyond the traditional endoreversible thermodynamics, which considers internal subsystems as reversible, i.e., without internal entropy production. Second, the mentioned work [35] goes beyond the use of paradigmatic, simple models and thus sheds light on the modeling characteristics of endoreversible systems in relation to real running heat engines.

Despite the extensive literature, the possibility of inducing optimization from stability in finite-time thermodynamics has not been studied. This is one of the goals of the present paper, which is focused on irreversible Carnot-like heat engine models with linear and inverse heat transfer laws, as two representative examples, assuming that the maximum power state behaves as a steady state. The kind of perturbations assumed in this analysis are from external sources or variations in the external control, such as those stemming from variations in the velocity of a piston. A benefit of this study is that, unlike the low-dissipation scheme, based on entropy production, the finite-time thermodynamics scheme allows a more straightforward analysis of its consequences on the working fluid and design of the heat engine, as it accounts for explicit heat transfer laws.

This paper is structured as follows: In Section 2, an overview of the endoreversible model is presented, along with some results on the maximum power regime for both the Newton and the phenomenological heat transfer laws. In Section 3, a multiobjective optimization is realized and the Pareto front is faced with endoreversible behavior. In Section 4, the stability dynamics are obtained from basic assumptions and the operation and relaxation times are compared in order to achieve a useful stability. In Section 5, the relaxation trajectories are analyzed in the space of the control variables as in the energetic space composed by power output, efficiency, and entropy generation. Finally some concluding remarks are presented.

## **2. A Quick Look at the Endoreversible Model**

The most basic endoreversible model consists on a baseline Carnot cycle whose working fluid operates with effective temperatures *Thw* and *Tcw* (isotherms of the Carnot cycle), irreversibly connected to two external reservoirs at temperature *Th* and *Tc* (*Th* > *Thw* > *Tcw* > *Tc*), as depicted in Figure 1.

**Figure 1.** Schematic representation of an endoreversible heat engine. The working fluid realizes a Carnot cycle operating between the isothermal processes at effective temperatures *Thw* and *Tcw* < *Thw*. The working fluid is irreversibly coupled to external reservoirs at temperatures *Th* and *Tc* < *Th*.

The endoreversible scheme requires knowledge of the heat transfer laws to model the heat fluxes between the external reservoirs and the working fluid. In a somewhat general case of heat transfer laws, *Qh* and *Qc* are expressed as:

$$Q\_{\rm h} = \left. \sigma\_{\rm hc} \sigma\_c \left( T\_h^k - T\_{\rm hw}^k \right) \operatorname{sgn} \left( k \right) t\_{\rm h} > 0, \tag{1}$$

$$Q\_{\mathbb{C}} = \left. \sigma\_{\mathbb{C}} \left( T\_{\mathbb{C}}^{k} - T\_{\text{cw}}^{k} \right) \text{sgn} \left( k \right) t\_{\mathbb{C}} < 0, \tag{2}$$

where *k* = 0 is the exponent of the heat transfer law (*k* = 1 refers to the known Newton law, *k* = −1 is known as the phenomenological law, *k* = 4 is the Stefan–Boltzmann law, and so on); *σ<sup>c</sup>* and *σ<sup>h</sup>* are the thermal conductances with units of J·K−*k*/s, here considered as constant since only small fluctuations around the operation regime are addressed. The function *sgn* (*k*) is defined as:

$$\text{sgn}\,(k) = \begin{cases} -1 & \text{if } \begin{array}{c} if \ k>0\\ -1 & \text{if } \begin{array}{c} k<0 \end{array} \end{array} \end{cases}$$

which is introduced for consistency with the convention that the heat flux entering the working fluid is positive and that the outgoing heat is negative; *th* and *tc* are the times at which the working fluid is in contact with the external heat reservoirs *Th* and *Tc*, respectively. In Equation (1), the *σhc* ≡ *σh*/*σ<sup>c</sup>* ratio is introduced, which is an important ingredient to obtain the upper and lower bounds of the efficiency of a given operation regime. The endoreversible hypothesis states that the entropy generation in the inner part of the engine is zero. Mathematically, this means that:

$$\frac{Q\_h}{T\_{hw}} + \frac{Q\_c}{T\_{cw}} = 0,\tag{3}$$

and from this constraint, the ratio of the operation time is given as:

$$\frac{t\_c}{t\_h} = \frac{T\_h^k - T\_{hw}^k}{T\_{cw}^k - \pi^k T\_h^k} \frac{T\_{cw}}{T\_{hw}} \sigma\_{hc} \tag{4}$$

where *τ* ≡ *Tc*/*Th*. With the endoreversible hypothesis constraining *tc*, in Equations (3) and (4), the efficiency *η*, power output *P*, and total entropy generation *S* become functions of {*Th*, *Tc*, *Thw*, *Tcw*, *σhc*, *th*, *σc*}. By using Equations (1) and (2), the explicit equations of these key thermodynamic magnitudes are the following:

$$\begin{aligned} \eta &= 1 + \frac{Q\_k}{Q\_h} = \begin{cases} 1 - \frac{T\_{\text{ctr}}}{T\_{\text{hur}}}, & k = \{1, -1\} \\ \frac{\sigma\_c \sigma\_{\text{hv}} (T\_{\text{ctr}} - T\_{\text{hv}}) (T\_h - T\_{\text{hv}}) (T\_{\text{ctr}} - \tau T\_h) \text{sgn}(k)}{T\_{\text{ctr}} T\_{\text{hv}} (\sigma\_{\text{hv}} - 1) + T\_h (\tau T\_{\text{hv}} - \tau \tau T\_h) \text{sgn}(k)} & k = 1 \end{cases} \\\\ \frac{\sigma\_c \sigma\_{\text{hv}} (T\_{\text{ctr}} - T\_{\text{hv}}) (T\_h - T\_{\text{hv}}) (T\_{\text{ctr}} - \tau T\_h) \text{sgn}(k)}{T\_h \left( T\_{\text{hv}}^2 (\tau\_{\text{ctr}} - \tau T\_h) + T\_{\text{hv}}^2 (T\_h - T\_{\text{hv}}) \tau \tau\_h \right)} & k = -1 \end{aligned} \tag{5}$$

$$\mathcal{S} = -\frac{Q\_h}{T\_h} - \frac{Q\_c}{T\_c} = \begin{cases} \frac{t\_h \sigma\_c \sigma\_{\text{hv}} (T\_h - T\_{\text{hv}}) (T\_{\text{ctr}} - \tau T\_{\text{hv}}) \text{sgn}(k)}{T\_h \hbar \tau \tau} & k = 1 \\\\ \frac{t\_h \sigma\_c \sigma\_{\text{hv}} (T\_h - T\_{\text{hv}}) (\tau T\_{\text{hv}} - \tau T\_{\text{ctr}}) \text{sgn}(k)}{T\_h^2 T\_{\text{hv}}^2 \tau} & k = -1 \end{cases} . \tag{6}$$

As can be seen in the definition of the power output, instantaneous adiabatic processes were considered. In order to deal with non-instantaneous adiabatics, a more detailed analysis of the working fluid and the geometry of the device needs to be made (for example, see [36]). In some works, these features have been analyzed by showing their influence in the power output and efficiency, but under certain circumstances, such as the large compression ratios in a piston, these times produce small effects on *η* and *P*, remaining compatible with the quasistatic nature of the internal process.

Notice that only *S* is proportional to *th*. As is shown later in Section 4 (see Equation (15)), the relaxation times are also proportional to this partial time; thus, a relationship between relaxation time and total entropy generation can be found. By fixing the value of *S*, a time scale for stability can be established. Besides this time scale (establishing the speed of the restitution dynamics), it is possible to analyze the energetic consequences of stability independently of *th* if one works with entropy production per cycle time, defined as *<sup>S</sup>*˙ <sup>≡</sup> *<sup>S</sup>*/(*tc* <sup>+</sup> *th*) = *<sup>S</sup>*/*th*(<sup>1</sup> <sup>+</sup> *tc*/*th*), taking advantage of the constraint on the total operation time through Equation (4).

The maximum power (MP) regime is specified by the temperatures *T*∗ *cw* resulting from the constraint *∂P*/*∂Tcw* = 0 and *T*<sup>∗</sup> *hw* under the condition *∂P*/*∂Thw* = 0. The resulting values are given by:

$$T\_{\rm hw}^{\*} \quad = \frac{T\_{\rm h}(\sqrt{\tau} + \sqrt{\sigma\_{\rm hc}})}{1 + \sqrt{\sigma\_{\rm hc}}} \quad k = 1,\tag{6}$$

$$T\_{\mathbb{C}\mathbb{D}}^{\*}\quad = \frac{T\_{h}\left(\tau+\sqrt{\sigma\_{h\varepsilon}}\right)}{1+\sqrt{\sigma\_{h\varepsilon}}}\quad k=1,\tag{7}$$

$$T\_{hw}^{\*} \quad = \frac{2T\_h \tau \left(1 + \sqrt{\sigma\_{hc}}\right)}{1 + \tau + 2\tau\sqrt{\sigma\_{hc}}} \quad k = -1,\tag{8}$$

$$T\_{\rm cw}^{\*} \quad = \frac{2T\_{\rm h}\tau \left(1 + \sqrt{\sigma\_{\rm h}}\right)}{2 + \sqrt{\sigma\_{\rm h}}(1 + \tau)} \quad k = -1. \tag{9}$$

From the optimal *T*∗ *cw* values, parametrization of *P*, *η*, and *S*˙ can be made (see dashed curve in Figure 2). The obtained parametric curve exhibits a parabolic-like behavior, where the values of *η* vary from 0 to *η<sup>C</sup>* ≡ 1 − *τ*, the Carnot efficiency (both limits corresponding to a zero power output), and with a single maximum in *P*. This is what is known in the literature as an endoreversible behavior, a kind of signature of the model. Some relevant features regarding the efficiency at maximum power from these endoreversible models with *k* = 1 and *k* = −1 are pointed out:


efficiency as in the low-dissipation model, where the self-optimization property has been studied. In this case (*k* = −1), *ηMP* is *σhc*-dependent, bounded by *ηMP* ∈ *η<sup>C</sup>* <sup>2</sup> , *<sup>η</sup><sup>C</sup>* 2−*η<sup>C</sup>* , according to whether *σhc* varies from 0 to ∞ [7].

**Figure 2.** Parabolic behavior of the *η*, *P*, and *S*˙ curve typical of the endoreversible model for (**a**) the *k* = 1 case and (**b**) the *k* = −1 case. The values *Th* = 500 K, *σhc* = 1, and *τ* = 0.4 weer fixed. Additionally, *σc* was chosen (for representation/comparison purposes) in such a way that the maximum power (MP) was 70 W in both cases, so that *η*, *P*, and *S*˙ ranged in similar intervals.

## **3. The Relevant Region for Optimization: The Pareto Front**

Before exploring stability dynamics, it is convenient to introduce a multiobjective optimization of the model. When looking for the best compromise among a variety of objective functions, the result is the so-called Pareto front (in the space of energetic functions) and the corresponding Pareto optimal set (in the space of operation variables).

To this end, a sorting algorithm is used by applying the concept of dominance [38]: A vector *v* = (*v*1, ... , *vn*) dominates another one *w* = (*w*1, ... , *wn*) if, and only if, *vi* ≥ *wi* ∀ *i* ∈ {1, ... , *n*} (if one is looking for a maximum, ≤ for a minimum) and there is at least one *j*, such that *vj* > *wj*, if one is interested in those vectors that are not dominated by any other. In other words, vectors that have the best value at least in one objective function. In this case, such a vector is formed by *η*, *P*, and *S*˙. The algorithm introduced here is a modification of the one introduced in [1–3], as follows:


In the present analysis, and in order to ensure convergence in the results, Kullback–Leibler (KL) divergence was introduced as a measure of the relative entropy *DKL* [39]. *DKL* was calculated between the probability distribution of the efficiencies of the Pareto front in the *i*th and the *i* − 1th iterations. As the probability distribution of *η* converges to the true Pareto front distribution, the entropy of the distribution converges as well. In such a case, the relative entropy *DKL* tends toward zero. The radius to extend the search region in the phase space decreases with the *DKL* value. When this relative entropy

is very small, there is no information gain in iterating more times; then, the search for new points in the Pareto optimal set stops. As mentioned before, *DKL* provides a measure of statistical convergence by indicating how distant two distributions are. If *DKL* = 0, then the information stemming from both distributions is the same. This is a relevant issue to demonstrate that the obtained trend is not due to the lack of additional iterations.

To obtain *DKL*, the range of possible values of *<sup>η</sup>*, each iteration was divided into <sup>√</sup>*<sup>N</sup>* (rounded to the upper next integer [40], with *N* being the number of random points added to the search in each iteration) equal intervals (or bins). In this way, the same partition was used to compute the discrete probability distributions, *ρ<sup>k</sup>* (for the *k* iteration). *DKL* was calculated by comparing *ρk*−<sup>1</sup> with *ρk*. *DKL*,*<sup>k</sup>* is given by:

$$D\_{KL,k}(\rho\_{k-1}||\rho\_k) = -\sum\_{i} \rho\_{k-1,i} \log \left(\frac{\rho\_{k,i}}{\rho\_{k-1,i}}\right),\tag{10}$$

allowing to determine how much information is gained by narrowing the search. In Figure 2, the Pareto front is depicted (green points), along with the endoreversible curve for *η*, *P*, and *S*˙ (dashed line), obtained from the parametric elimination of *Thw* (*Tcw* is constant for the fixed values of the *σhc*, *σc*, *τ*, and *Th* parameters; see Equations (7) and (9)). Notice that even when the endoreversible curve depends on a first constraint (*Tcw* satisfies *∂P*/*∂Tcw* = 0), the Pareto front has nothing to do with. Nonetheless, the Pareto front lies over the endoreversible curve, covering the regions from maximum power to maximum efficiency and minimum entropy production. This is consistent with the literature, where the mentioned part of the curve has been denoted as the relevant region for optimization. For optimization purposes, additional figures of merit as compromise functions between these three (such as the Ecological function [41] or the Omega function [42]) will not provide further information to the Pareto front. After all of these considerations, analysis of stability dynamics can be made.

## **4. Stability Dynamics and Relaxation Times**

One goal of this analysis was to recover the previous results obtained in the low-dissipation scheme (with no direct connection to working fluid particularities) and to lay the groundwork for a more direct connection with the properties of the working fluid and design parameters. The perturbations we were interested in involved only those of an external nature. The temperature of the external reservoirs remained constant, as well as the conductances of the heat transfer laws (lastly depending on temperature) for only small variations in the operation regime. To illustrate these points, consider a 2-D piston describing a Carnot cycle. The velocity of the piston will determine the effective temperatures of the isotherms of the particles inside the piston (see [32]). The external control (and its energetic cost) needs to remain constant as the velocity of the piston is not accounted for in the thermodynamic description of the gas inside of it. However, fluctuations in this velocity, as well as possible cyclic variability, will lead to variations in the effective temperatures *Thw* and *Tcw*. Then, a compromise between the energetic cost of the control and the thermodynamic consequences due to stability is of interest.

A simple approach to tackle this problem is the following: Since the operation regime is entirely defined by *Thw* and *Tcw*, the dynamics involving only these two variables are to be addressed. Fluctuations in these temperatures will affect the heat exchanged between the working fluid and the external reservoirs, *Qh* and *Qc*. From a Taylor expansion in the first order, *Qh* and *Qc* near the MP state can be written in matrix form as:

$$
\begin{pmatrix} Q\_h - Q\_h^\* \\ Q\_c - Q\_c^\* \end{pmatrix} = \begin{pmatrix} \frac{\partial Q\_h}{\partial T\_{hw}}\Big|\_\ast & \frac{\partial Q\_h}{\partial T\_{cw}}\Big|\_\ast \\\ \frac{\partial Q\_c}{\partial T\_{hw}}\Big|\_\ast & \frac{\partial Q\_c}{\partial T\_{cw}}\Big|\_\ast \end{pmatrix} \cdot \begin{pmatrix} T\_{hw} - T\_{hw}^\* \\\ T\_{cw} - T\_{cw}^\* \end{pmatrix} = \begin{pmatrix} a & 0 \\ \ \beta & \gamma \end{pmatrix} \cdot \begin{pmatrix} T\_{hw} - T\_{hw}^\* \\\ T\_{cw} - T\_{cw}^\* \end{pmatrix} \tag{11}
$$

$$
\Rightarrow \quad \begin{pmatrix} T\_{hw} - T\_{hw}^\* \\ T\_{cw} - T\_{cw}^\* \end{pmatrix} = - \begin{pmatrix} \frac{1}{\bar{a}} & 0 \\ -\frac{\beta}{a\gamma} & \frac{1}{\bar{\gamma}} \end{pmatrix} \cdot \begin{pmatrix} Q\_h^\* - Q\_h \\ Q\_c^\* - Q\_c \end{pmatrix}, \tag{12}
$$

where *α*, *β*, and *γ* are the elements of the Jacobian matrix evaluated in the MP state (denoted by ∗). It is noted again that this analysis assumes that the MP state is a steady stationary state. Within the first order scheme [43], the simplest relationship for an autonomous system for the two variables *Thw* and *Tcw* is:

$$
\begin{pmatrix} \dot{T}\_{hw} \\ \dot{T}\_{cw} \end{pmatrix} = - \begin{pmatrix} A & 0 \\ 0 & B \end{pmatrix} \cdot \begin{pmatrix} T\_{hw} - T\_{hw}^\* \\ T\_{cw} - T\_{cw}^\* \end{pmatrix} \tag{13}
$$

which is a good approximation near a stable point. The coefficients *A* and *B* determine the restitution strength (with units of s−1). The bigger they are, the faster the system will return to the steady state. By substituting Equation (12) into Equation (13), with the resulting dynamics:

$$
\begin{pmatrix} \mathcal{T}\_{hw} \\ \mathcal{T}\_{cw} \end{pmatrix} = \begin{pmatrix} \frac{\mathcal{A}}{\alpha} & 0 \\ -\frac{B\mathcal{S}}{a\gamma} & \frac{\mathcal{B}}{\gamma} \end{pmatrix} \cdot \begin{pmatrix} Q\_h^\* - Q\_h \\ Q\_c^\* - Q\_c \end{pmatrix} . \tag{14}
$$

The magnitudes of the relaxation times *t*<sup>1</sup> = −1/*λ*<sup>1</sup> and *t*<sup>2</sup> = −1/*λ*<sup>2</sup> are obtained from the eigenvalues, *λ*1,2, of the square matrix in the right-hand side of Equation (14). For the cases *k* = 1 and *k* = −1, they are:

$$\begin{split} t\_1 &= \begin{cases} \frac{t\_k \sigma\_\ell \sigma\_{\mathcal{U}}}{A} & k = 1 \\\\ \frac{t\_k \sigma\_\ell \sigma\_{\mathcal{U}} \left(1 + \tau + 2\tau \sqrt{\sigma\_{\mathcal{U}}}\right)}{4kT\_h^2 \tau^2 \left(1 + \sqrt{\sigma\_{\mathcal{U}}}\right)^2} & k = -1, \end{cases} \\ t\_2 &= \begin{cases} \frac{t\_k \sigma\_\ell \sigma\_{\mathcal{U}} \left(1 - \sqrt{\tau}\right)}{B\left(\sqrt{\sigma\_{\mathcal{U}}} + \sqrt{\tau}\right)} & k = 1 \\\\ \frac{t\_k \sigma\_\ell \sigma\_{\mathcal{U}} \left(1 + \tau + 2\tau \sqrt{\sigma\_{\mathcal{U}}}\right) \left(1 - \tau\right)}{4kT\_h^2 \tau^2 \left(1 + \sqrt{\tau\_{\mathcal{U}}}\right)^2} & k = -1. \end{cases} \end{split} \tag{15}$$

Due to the units of the two matrices, a K/J unit factor should be accounted for in the final expression on the right-hand side of Equations (15), providing the correct units for the relaxation times (in s). As mentioned before, *A* and *B* provide the strength of the restitution dynamics, i.e., if they have large values, then the relaxation times are short. Notice that both the relaxation times and the total entropy generation are proportional to the contact time *th*; then, *S* can be used to define a characteristic time scale for relaxation. Since there are no reasons to provide a preferred relaxation in the heat exchange *Qc* or *Qh*, it can be assumed that *t*<sup>1</sup> = *t*<sup>2</sup> (this requirement can be tuned according to the specific conditions of a heat device at hand). Additionally, when looking for stability of a cyclic process, it is desirable that the stability is achieved in times shorter than those of the operation time; that is:

$$t\_1 + t\_2 = 2t\_1 \le t\_h + t\_c = t\_h \left( 1 + \frac{t\_c}{t\_h} \right). \tag{16}$$

From this constraint and Equation (4), the relaxation times are bounded by:

$$t\_1 = t\_2 \le \begin{cases} \frac{t\_h}{2} \left( 1 + \sqrt{\sigma\_{hc}} \right) & k = 1 \\\\ \frac{t\_h \left( 1 + \sqrt{\sigma\_{hc}} \right) \left( 1 + \tau \sqrt{\sigma\_{hc}} \right)}{2 + \sqrt{\sigma\_h \epsilon} \left( 1 + \tau \right)} & k = -1 \end{cases} \tag{17}$$

where the equality corresponds to the case where *tc* + *th* = *t*<sup>1</sup> + *t*2. From now on, equality is assumed. Then, the resulting dynamics (Equation (14)) for *k* = 1 are:

$$\begin{array}{rcl} \dot{T}\_{\text{lfw}} &=& -\frac{2r\_{l}\sigma\_{\text{lbf}}\left(T\_{\text{lw}}\left(1+\sqrt{\sigma\_{\text{lw}}}\right) - T\_{\text{l}}\left(\sqrt{\sigma\_{\text{lw}}}+\sqrt{\tau}\right)\right)}{\left(1+\sqrt{\sigma\_{\text{lw}}}\right)^{2}},\\ \dot{T}\_{\text{cw}} &=& -\frac{2r\_{c}\sigma\_{\text{lw}}T\_{\text{l}}\sqrt{\tau}}{1+\sqrt{\sigma\_{\text{lw}}}} \left(\frac{T\_{\text{lw}}\left(1+\sqrt{\sigma\_{\text{lw}}}\right)}{T\_{\text{l}}\left(\sqrt{\sigma\_{\text{lw}}+\sqrt{\tau}}\right)} + \frac{T\_{\text{lw}}}{T\_{\text{l}}\sqrt{\tau}}\left(\frac{T\_{\text{h}}}{T\_{\text{lw}}}-1\right) - \frac{1-\tau+2\left(\sqrt{\sigma\_{\text{lw}}}+\sqrt{\tau}\right)}{\left(1+\sqrt{\sigma\_{\text{lw}}}\right)\left(1+\sqrt{\tau}\right)}\right), \end{array} \tag{18}$$

and for *k* = −1,

$$\begin{array}{rcl} \dot{T}\_{hw} &=& -\frac{4\sigma\_{\mathcal{C}}\sigma\_{\mathcal{U}}}{T\_{hw}} \cdot \frac{\left(\frac{T\_{\text{flux}}}{T\_{\text{k}}}\left(\frac{1\pm\varepsilon}{2\overline{\tau}} + \sqrt{\sigma\_{\text{k}}}\right) - 1 - \sqrt{\sigma\_{\text{k}}}\right)\left(1 + \sqrt{\sigma\_{\text{k}}}\left(\frac{1\pm\varepsilon}{2}\right)\right)}{\left(1 + \sqrt{\sigma\_{\text{k}}}\right)^{2}\left(1 + \tau\sqrt{\sigma\_{\text{k}}}\right)} \\\\ \dot{T}\_{\text{cW}} &=& -\frac{4\sigma\_{\mathcal{C}}\sigma\_{\text{k}}}{T\_{hw}} \cdot \frac{\frac{\tau T\_{\text{flux}}}{I\_{h}}\left(\frac{1\pm\varepsilon}{2\overline{\tau}} + \sqrt{\sigma\_{\text{k}}}\right)^{2} - \left(1 + \sqrt{\sigma\_{\text{k}}}\right)\left(1 + \tau\sqrt{\sigma\_{\text{k}}} + \frac{T\_{\text{aux}}}{I\_{h}}\left(1 + \sqrt{\sigma\_{\text{k}}}\left(\frac{1\pm\varepsilon}{2}\right)\right)\right) + \frac{T\_{\text{aux}}}{I\_{hw}}\left(1 + \sqrt{\sigma\_{\text{k}}}\right)\left(1 + \tau\sqrt{\sigma\_{\text{k}}}\right)}. \end{array} \tag{19}$$

From these dynamic equations, it is possible to calculate a relaxation velocity, *v* = *T*˙ 2 *hw* <sup>+</sup> *<sup>T</sup>*˙ <sup>2</sup> *cw* (in K/s), which indicates how fast the system is evolving toward the steady state. In Figure 3, iso-velocity contours are depicted for the *k* = {1, −1} cases. Two trajectories are represented to provide an idea of how fast the system can return to the MP state. In both cases, the operation time is indicated, along with the time of evolution of each curve. As can be seen, 1 or 2 s are enough to drive the system close to the stationary state; meanwhile, the cycle time is approximately 300 s. Another feature in Figure 3 is that far from the steady state, the system evolves faster. Different marks at equal time intervals are placed over the trajectories, and in every case, the system travels a longer distance in the first interval. In both cases, the *Tcw* direction is slower. It can be confirmed that the depicted trajectories evolve more quickly in the horizontal direction, while the final approach is mostly in the vertical direction. Notice also that the *k* = −1 case has the fastest dynamics.

**Figure 3.** Isocontours of the relaxation velocity for (**a**) the *k* = 1 case and (**b**) the *k* = −1 case. The values *Th* = 500 K, *σhc* = 1, and *τ* = 0.4 were fixed. As in Figure 2, *σ<sup>c</sup>* and *th* were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases. By fixing *S* at MP conditions, a time scale for *th* was established, and therefore, a scale for the relaxation time. The qualitative behavior for the other values of *S* was similar to the one presented here.

## **5. Thermodynamics of the Relaxation Trajectories**

From the analysis of the evolution of trajectories toward the steady state, it is not obvious which energetic implications arise. To address this issue, the dynamic equations were numerically solved. Several trajectories with a starting point within a radius of 30 K (in the *Thw*–*Tcw* space) from the MP state were analyzed. Then, these trajectories were mapped into the energetic space (*η*, *P*, and *S*˙).

As it can be checked in Equations (4), (17) and (18), the conductance ratio *σhc* plays a key role both in the dynamic equations and in the energetic properties of the system. In this regard, it is noted that the limits *σhc* → {0, ∞} are not physical from a dynamic point of view. For both the *k* = {−1, 1} cases, the limit *σhc* → 0 leads to null power output *P*, entropy *S* (Equation (5)), and also null restitution strengths (see Equation (18)), i.e., {*P*, *S*, *A*, *B*} → 0; thus, there is no power output and there are no restitution forces. In comparison, the limit *σhc* → ∞ leads to {*S*, *A*, *B*} → ∞; thus, there is an infinite entropy generation and infinite restitution forces. For this reason, in the present analysis, the more realistic cases *<sup>σ</sup>hc* <sup>=</sup> {10−6, 1, 106} were considered.

In Figure 4, the four rows display the trajectories obtained from numerically solving the dynamics in Equation (18) for the *k* = 1 case. The first row displays the trajectories in the *Thw* and *Tcw* space, while the next rows show, respectively, the mapping of the same trajectories over the *η*–*P*, *S*˙–*P*, and *S*˙–*η* spaces. The three columns are for the abovementioned representative conductance values: (a) *σhc* = 10−6, (b) *σhc* = 1, and (c) *σhc* = 106. In the first row, according to the initial state in each trajectory, a clock hour analogy was used; in this way, the purple line corresponds to 12 o'clock and red line to 3 o'clock. The Pareto optimal set is displayed (green points), together with some iso-velocity contours. The white curve indicates the position of the MP state as *σhc* varies from 0 to ∞. All of these points produce the same efficiency at MP given by *η*<sup>∗</sup> = *ηCAN*. For the rest of the rows, the Pareto front and the endoreversible curve are presented as well. Note also that in the three configurations (columns), the total operation time is the same, but the *tc*/*th* (and *σc*) ratio varies, as can be seen in the legend in each case. This figure reveals several key features:

(a) For a small *σhc* (left column), only trajectories between 9 and 3 o'clock are present in the *Thc*–*Tcw* space (denoted in colors ranging from red to blue). Most of the observed trajectories (especially those with darker colors) evolve to the stable state near the endoreversible curve in such a way that the power and efficiency progressively increase their values, while entropy production decreases (see rows 2–4 in the first column). In other words, the relaxation process mainly drives the system toward a thermodynamic steady state, thereby enhancing the thermodynamic performance of the engine (*η* and *P* increase and *S*˙ decreases).

(b) For *σhc* = 1 (second column), the trajectories arriving at the stable point from all directions are clearly observed. Those in darker colors (from 9 to 3 o'clock) evolve, as in the above case, but closer to the endoreversible curve, showing an increase in power and efficiency and a decrease in entropy generation. In contrast, the trajectories denoted in colors ranging from orange to green (between 3 and 7 o'clock) evolve to the steady state near the Pareto front, with an increase in power output but a decrease in efficiency and higher entropy production. Note also in this configuration how the relaxation is well balanced between the trajectories approaching the endoreversible limit on the Pareto front side and to the other side. However, the curves near the Pareto front are longer, meaning that random perturbations tend to favor a locus directed toward the Pareto front.

(c) For a large *σhc* (last column), the trajectories evolve directly toward the endoreversible curve first. In this part of the trajectory, *η*, *P*, and *S*˙ are enhanced simultaneously. Later on, the system evolves to the steady state, either through the endoreversible limit or the Pareto front.

(d) Although the sizes of the perturbations in the *Thw*–*Tcw* space are the same, the case with the small *σhc* allows larger variations of power output, while the case where *σhc* = 1 exhibits the smallest fluctuations in the *η*–*S*˙ plane.

(e) The relaxation times (Equation (17)) are directly proportional to *th*, and from the expressions of *P*, *η*, and *S* (see Equation (5)), only entropy generation depends on *th*. Thus, the characteristic time scale of the relaxation is linked to the entropy scale of the system, a feature that connects the stability with the thermodynamics of the system. In this entropy–control point of view, the dynamics for the small values of *σhc* in the left column apply to greater *th* values and, as a consequence, when the contact time of the heat engine with the cold reservoir is small. On the contrary, the dynamics for the large values of *σhc* in the right column apply to smaller *th* values, i.e., when the contact time of the heat engine with the cold reservoir is large.

**Figure 4.** Trajectories in the *Thw*–*Tcw* space and the mapping over the *η*–*P*, *S*˙–*P*, and *S*˙–*η* spaces for *k* = 1. (**a**) the *σhc* = 10−<sup>6</sup> case; (**b**) the *σhc* = 1 case; (**c**) the *σhc* = 10<sup>6</sup> case. The values *Th* = 500 K and *τ* = 0.4 were fixed. As in Figure 2, *σ<sup>c</sup>* and *th* were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases.

In Figure 5, the corresponding configurations for *k* = −1 are depicted. Notice that the behavior of the trajectories is almost the same as that in the *k* = 1 case, but here, the maximum power stationary state is linked to different efficiencies as *σhc* increases from 0 to ∞. In particular, when *σhc* → 0, the efficiency at maximum power is *η*<sup>∗</sup> = <sup>1</sup>−*<sup>τ</sup>* <sup>2</sup> <sup>=</sup> *<sup>η</sup><sup>c</sup>* <sup>2</sup> ; when *<sup>σ</sup>hc* <sup>→</sup> 1, *<sup>η</sup>*<sup>∗</sup> <sup>=</sup> <sup>2</sup>(1−*τ*) <sup>3</sup>+*<sup>τ</sup>* <sup>=</sup> <sup>2</sup>*η<sup>C</sup>* <sup>4</sup>−*η<sup>C</sup>* ; when *<sup>σ</sup>hc* <sup>→</sup> <sup>∞</sup>, *<sup>η</sup>*<sup>∗</sup> <sup>=</sup> <sup>1</sup>−*<sup>τ</sup>* <sup>1</sup>+*<sup>τ</sup>* <sup>=</sup> *<sup>η</sup><sup>C</sup>* <sup>2</sup>−*η<sup>C</sup>* . The dynamics for different *σhc* are no longer linked to equal operation times *tc* + *th*. Another difference is that the relaxation times are noticeably smaller and, therefore, the restitution forces are stronger. Additionally, as *σhc* decreases (*η*<sup>∗</sup> also decreases), the same temperature perturbations lead to the farthest starting points in the energetic planes. In this way, the case of tlarge *σhc* (with a high efficiency) is the one with the smallest

drops in *η*, *P*, and *S*˙ (with shorter trajectories). This suggests that the system becomes more stable as the efficiency in the steady state increases.

**Figure 5.** Trajectories in the *Thw*–*Tcw* space and mapping over the *<sup>η</sup>*–*P*, *<sup>S</sup>*˙–*P*, and *<sup>S</sup>*˙–*<sup>η</sup>* spaces for *<sup>k</sup>* <sup>=</sup> <sup>−</sup>1. In column (**a**) the *σhc* = 10−<sup>6</sup> case; (**b**) the *σhc* = 1 case; (**c**) the *σhc* = 10<sup>6</sup> case. The values *Th* = 500 K and *τ* = 0.4 were fixed. As in Figure 2, *σ<sup>c</sup>* and *th* were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases.

## **6. Concluding Remarks**

It was shown that the Pareto front, which represents the best compromise among power output, efficiency, and entropy generation is related to endoreversible behavior, obtained herein by analyzing the Newton and phenomenological heat transfer laws in the context of Carnot-like models. These findings corroborate those obtained in the low-dissipation scheme.

A set of dynamic equations were found based only on the assumption that the endoreversible heat engine has a stationary state and from a Taylor expansion of the input/output heat fluxes in the first order. A relationship between relaxation times and total operation time was obtained. It was shown that the trajectories that lead the system back to the stationary state require much shorter times than the cycle time, allowing the system to work under continuous cyclic operation.

Mapping of the relaxation trajectories into the energetic space allowed for an analysis of the performance consequences of stability. The degree of symmetry of the conductance ratio in the input/output heat exchange is the main valuable factor for stability dynamics. For small *σhc* values, the thermodynamic trajectories improve the efficiency and power values with a decrease in entropy generation, evolving near the endoreversible behavior. For larger *σhc* values, the first part of the thermodynamic trajectories improve the efficiency and power values with a decrease in entropy generation, evolving toward the endoreversible behavior and the Pareto front. In between these two situations, i.e., for equal conductance values, (*σhc* = 1) trajectories with decreasing efficiencies and increasing entropy generation can be found, with a preference for evolving near the Pareto front.

In summation, the stability of the irreversible Carnot-like heat engine exhibits two interesting behaviors in which the fluctuations around the stationary state (due to external perturbations) would likely maintain the system in an optimum state, or produce self-optimization induced by the stability. A biased control of the operation parameters could result in an economic saving, as energy is needed for the fine-tuning of parameter control. Finally, in this work the thermodynamic functions that were selected as the most relevant for this heat engine model were *η*, *P*, and *S*˙. It is very likely that by analyzing different stationary states, such as those predicted by the compromise Ecological [41] or Omega [42] operation performances, the role of the Pareto front would be more evident, as occurred in the case of the low-dissipation heat engine and refrigerator [1,2].

Another valuable remark is that the total entropy generation defines a time scale for both the operation and relaxation times. Finally, this analysis laid down the grounds to analyze heat devices dependent on working fluid properties, where the endoreversible hypothesis plays a relevant role. Especially when the thermalization mechanisms are fast enough compared to the cycle time in agreement with local equilibrium. This allows to incorporate, in a straightforward way, stochastic-type perturbations into the analysis through additive noise in Equation (14).

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

**Funding:** Junta de Castilla y León, Project No. SA017P17; J.G.A. acknowledges partial financial support from Instituto Universitario de Física Fundamental y Matemáticas (IUFFyM) at Universidad de Salamanca.

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

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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/).
