Optimization-based investigations of a two-phase thermofluidic oscillator for low-grade heat conversion

Abstract

Background

The non-inertive-feedback thermofluidic engine (NIFTE) is a two-phase thermofluidic oscillator capable of utilizing heat supplied at a steady temperature to induce persistent thermal-fluid oscillations. The NIFTE is appealing in its simplicity and ability to operate across small temperature differences, reported as low as 30 ^∘C in early prototypes. But it is also expected that the NIFTE will exhibit low efficiencies relative to conventional heat recovery technologies that target higher-grade heat conversion. Mathematical modeling can help assess the full potential of the NIFTE technology.

Results

Our analysis is based on a nonlinear model of the NIFTE, which we extend to encompass irreversible thermal losses. Both models predict that a NIFTE may exhibit multiple cyclic steady states (CSS) for certain design configurations, either stable or unstable, a behavior that had never been hypothesized. A parametric analysis of the main design parameters of the NIFTE is then performed for both models. The results confirm that failure to include the irreversible thermal losses in the NIFTE model can grossly overpredict its performance, especially over extended parameter domains. Lastly, we use the NIFTE model with irreversible thermal losses to assess the optimization potential of this technology by conducting a multi-objective optimization. Our results reveal that most of the optimization potential is achievable via targeted modifications of three design parameters only. The Pareto frontier between exergetic efficiency and power output is also found to be highly sensitive to these optimized parameters.

Conclusions

The NIFTE is of practical relevance to a range of applications, including the development of solar-driven pumps to support small-holder irrigation in the developing world. Given its low capital cost, potential improvements greater than 50% in efficiency or power output are significant for the uptake of this technology. Conventional heat recovery technologies are known to have higher efficiencies than those reported in this work, but they also have more complex designs and operations, higher capital costs, and may not even be feasible for the temperature differences considered herein. Future work should focus on confirming this model-based assessment via dedicated experimental campaigns and on investigating design modifications to mitigate irreversible thermal losses.

Background

High-temperature heat engines are ubiquitous in the power generation industry, stationary distributed power (co)generation, the transport and many other sectors. This demand for high-grade heat has led to ever increasing rates of fossil-fuel consumption and resource depletion, while at the same time concerns exist in relation to the emissions associated with the combustion of these fuels. There is therefore a pressing need for exploiting alternative heat sources.

Broadly speaking, heat available at a relatively low temperature and for which there is no direct use is classified as “low-grade”. Such heat is abundant in the process and chemical industries, where it is often “wasted” by release into the environment without further utilization. Large amounts are also generated and available in solar collectors and geothermal wells, among other applications. Unfortunately, all technologies for low-grade heat conversion to power are associated with low thermal efficiencies [1, 2], a restriction imposed by the second law of thermodynamics [3].

Thermofluidic oscillators are a class of engines for converting heat into useful work via the creation of oscillatory fluid motion. The non-inertive-feedback thermofluidic engine (NIFTE) is a particular type of two-phase thermofluidic oscillator [4–6]. The description of the configuration and step-by-step operation of the NIFTE are given as part of Fig. 1. The working fluid exists in both the liquid and vapor phases, and changes phase periodically during the operation (evaporation and condensation). It can be tailored to suit the available external heat source and sink temperatures. A NIFTE is simple to construct and has few moving parts in comparison with conventional engines, making it both a reliable and low capital-cost technology, that is capable of converting low-grade heat.

Fig. 1

NIFTE fluid pump schematic, adapted from Reference [25]. Configuration. A NIFTE consists of four major elements: the displacer cylinder ➇ that contains the hot ➃ and cold ➄ heat-exchangers; the adiabatic vapor region ➈; the power cylinder ➆ that is connected to the load ➁; and the feedback connection tube ➅. The white area denotes the vapor phase, while the grey area denotes the liquid phase of the working fluid. Operational cycle. Low-grade heat from an external source is supplied to the system via the hot heat-exchanger (HHX) ➃ to vaporize the working fluid. The pressure in the adiabatic vapor region ➈ thus increases, causing the liquid level ➀ in the power cylinder ➆ to fall below the equilibrium level. The liquid level ➂ in the displacer cylinder ➇ drops as a result of the hydrostatic pressure via the feedback connection ➅. The vapor comes into contact with the cold heat-exchanger (CHX) ➄, where heat is removed from the system, causing the vapor to condense and the liquid height to rise back above the equilibrium position. Sustained oscillations in the load ➁ are thereby created by the cyclic condensation and evaporation of the working fluid in the cylinders

It has been established experimentally on early prototypes [7] that a NIFTE is capable of operating with a temperature difference as small as 30 ^∘C between the heat source and sink, and more recent versions of the technology may be able to operate with temperature differences even smaller than this. An early prototype was reported as achieving thermal efficiencies up to 1% and exergetic efficiencies—the thermal efficiency as a fraction of the Carnot efficiency—between 1–2%, when operated with a heat source at 65–90 ^∘C, a heat sink at 4–12 ^∘C, and n-pentane as the working fluid [6, 8].

Such efficiencies can be considered low in comparison to those of other heat conversion technologies that are typically designed and intended for higher-grade heat sources, such as organic Rankine cycle (ORC) systems, Stirling and thermoacoustic engines [3]. In general, thermal efficiencies in the approximate range of 10–20% are expected from ORC systems, also depending on the system size, the lower end achievable at heat source temperatures of about 100–150 ^∘C and the upper end at up to 400–500 ^∘C [9]; while Stirling and thermoacoustic engines have much lower efficiencies at those temperatures, but can reach thermal efficiencies between 20–30% at heat-source temperatures of 600–800 ^∘C [3]. In one example [10], a small-scale regenerative ORC power system comprising a scroll expander and with R123 as the working fluid was reported as having a thermal efficiency upwards of 7% at a heat-source temperature of 165–183 ^∘C and an ambient air temperature of 22.5 ^∘C cooling the condenser. In a traveling-wave thermoacoustic engine with air as the working gas [11], a thermal efficiency of 14.5% was reported under 25 bar, with a heat-source temperature of 650 ^∘C and a cooling temperature of 12 ^∘C. Importantly, it can be seen that these conventional power generation systems are neither technically suitable nor economically feasible when the heat-source temperature drops below 100 ^∘C. Furthermore, the construction and operation of these systems are invariably more complicated than those for a thermofluidic oscillator, which has no or few moving parts or mechanical components and dynamic seals.

An alternative type of two-phase thermofluidic oscillator, known as the Up-THERM engine, has also been recently investigated for heat conversion applications [12, 13]. A lab-scale Up-THERM engine had a predicted exergetic efficiency of about 0.8%, when operating with a heat source and sink at 450 ^∘C and 10 ^∘C, respectively, and n-pentane as the working fluid [14]. It is also worth noting that an optimized pilot-scale Up-THERM engine with the same working fluid was predicted to yield promising exergetic efficiencies between 10–40% for a heat source temperature in the range 200–400 ^∘C [13], which is close to the efficiency of its ORC counterparts and could even be higher with optimally selected working fluids [15].

Nearly all of the mathematical models proposed so far to analyze the dynamical behavior of the NIFTE feature linear dynamics. Smith [4–6] first developed a linear model, based upon earlier work on thermoacoustic engines [16] and electrical analogies for thermal and fluid processes [17, 18]. In this early work, the NIFTE components were described by linear equations in analogy with electrical components—namely, resistors to account for fluid drag and thermal resistance, and capacitors to represent hydrostatic pressure and vapor compressibility—and the heat-exchange rate between the hot and cold elements and the working fluid was assumed to be either linear or constant. A complete engine model then connects these electrical components into a circuit according to the actual thermal-fluid technology layout. Inductors were later introduced in the electrical analogy to account for fluid inertia [7, 19], which allowed improvements to the accuracy of the predictions, and an effort to consider performance over wide parameter ranges was also made [20].

Notwithstanding their ability to provide quantitative predictions of the NIFTE performance via numerical simulation at marginal stability, these linear models cannot describe the sustained and robust periodic oscillations of the working fluid as observed in practice, which in essence is a nonlinear phenomenon. Subject to a steady temperature difference between the two heat-exchangers, the engine reaches a cyclic steady state (CSS) after a sufficiently long settling time, where the amplitude and frequency of the oscillations remain constant. When variations in the heat-exchanger temperature difference or other external disturbances occur, which is inevitable in practice, the engine will eventually again return to a CSS after some time. This is the signature of an asymptotically stable limit cycle in the dynamics, which may only occur in nonlinear dynamic systems.

Inspired by earlier work on exploiting the nonlinearity of core physical elements to predict the efficiency and robustness of thermodynamic engines [21–24], Markides and co-workers [25] were the first to propose a nonlinear inertive model of the NIFTE, referred to as the NTP model. The key nonlinearity in this model corresponds to the physical saturation in heat transfer rate between the hot and cold elements and the working fluid, whereas the rest of the engine components are described using the same linear electrical analogy used in earlier work.

Failure to account for various irreversible thermal losses in the early NIFTE models can lead to an overestimation of the engine’s efficiency and work output by an order of magnitude or more. Smith [6] identified the two predominant exergy losses to be: (i) shuttle losses, associated with the irreversible two-way heat transfer between the working fluid and the cylinder walls; and (ii) entrance and retrograde condensation losses, caused by the alternating parasitic condensation and evaporation of the working fluid. Solanki and Markides [8, 26] were the first to account explicitly for these losses in a linear NIFTE model, which led to much more reliable efficiency predictions. However, a similar extension of the NTP model to encompass exergy losses has not yet been reported in the literature.

This paper describes a detailed, model-based assessment of the performance of the NIFTE technology. We base our assessment on the existing NTP model, and propose a model extension, referred to as NTP-TL, that encompasses the irreversible thermal losses. Our methodology entails a model-based analysis of key design parameters—via the resistors, capacitors, inductors and other gain values in the electrical analogy—in terms of exergetic efficiency and average power output under CSS operation, that leads to new insight into the operation and limits of the NIFTE. We also perform a multi-objective optimization based on the NTP-TL model to explore the optimal trade-off between efficiency and power output, and gain further insight into NIFTE’s optimization potential.

Mathematical modeling

The two main thermal-fluid transport processes that govern the operation of a NIFTE device are: (i) pressure differences, driving working-fluid flows through fluid components; and (ii) temperature differences, responsible for heat flows in thermal components. We provide a summary of the models developed to describe these processes in the fluid and thermal domains first, followed by a complete statement of both the NTP and NTP-TL device models, and a description of the performance indicators considered for the assessment. A nomenclature of the notation used for mathematical modeling can be found in Table 1. For further details about the model development, refer to the original paper [25].

Table 1 Nomenclature

Fluid domain

The main components ➁, ➅, ➆, ➇, ➈ of the engine (cf. labels in Fig. 1) are modeled based on electrical analogies. The drag experienced by the working fluid in the load ➁ (due to viscosity) and in the feedback connection ➅ (due to small diameter and/or a valve) are represented by resistors, R_i. In particular, the model assumes that the fluid flow in these elements is akin to viscous laminar flow in a smooth pipe. The inertive effect experienced by the oscillating liquid flow in the load ➁, feedback connection ➅, power cylinder ➆, and displacer cylinder ➇ are modeled by inductors, L_i. The hydrostatic pressures at the bottom of the power cylinder ➆ and displacer cylinder ➇, and the adiabatic compression/expansion in the vapor region ➈ are represented by capacitors, C_i. The processes of vapor compression and expansion undergone by the working fluid in the vapor region ➈ are assumed to be isentropic (i.e. both adiabatic and reversible).

These R,L,C parameters may be expressed in terms of the physical and geometric properties of the NIFTE, namely density, ρ and dynamic viscosity, μ of the working and pumped fluids (indexed with _wf,_pf, respectively), and diameter, d, length, l and cross-sectional area, A of the tubes and valves (indexed with _l,_f,_p,_d,_ad in reference to the components ➁, ➅, ➆, ➇, ➈, respectively):

$$\begin{array}{*{20}l} R_{\mathrm{l}} = & \frac{128\mu_{\text{pf}} l_{\mathrm{l}}}{\pi d_{\mathrm{l}}^{4}}\:, & \!\!\!\! R_{\mathrm{f}} = & \frac{128\mu_{\text{wf}} l_{\mathrm{f}}}{\pi d_{\mathrm{f}}^{4}}\:,\\ L_{\mathrm{l}} = & \frac{\rho_{\text{pf}} l_{\mathrm{l}}}{A_{\mathrm{l}}}\:, &\!\!\!\! L_{\mathrm{f}} = & \frac{\rho_{\text{wf}} l_{\mathrm{f}}}{A_{\mathrm{f}}}\:, &\!\!\!\! L_{\mathrm{p}} = & \frac{\rho_{\text{wf}} l_{\mathrm{p}}}{A_{\mathrm{p}}}\:, &\!\!\! L_{\mathrm{d}} = & \frac{\rho_{\text{wf}} l_{\mathrm{d}}}{A_{\mathrm{d}}},\\ C_{\mathrm{p}} = & \frac{A_{\mathrm{p}}}{\rho_{\text{wf}} g}\:, &\!\!\!\! C_{\mathrm{d}} = & \frac{A_{\mathrm{d}}}{\rho_{\text{wf}} g}\:, &\!\!\!\! C_{\text{ad}} = & \frac{V^{\circ}}{\gamma P^{\circ}}\:, \end{array} $$

(1)

with V^∘ and P^∘, the equilibrium (time-averaged) volume and pressure in the gas phase; γ, the adiabatic index; and g, the gravitational acceleration. All of the physical properties, and therefore the R,L,C parameters themselves, are assumed to be temperature independent. But for a particular fluid, the R,L,C parameter values can be adjusted by varying the geometric properties (length l, diameter d, and/or cross-sectional area A) of the corresponding elements, which provides the basis for the parametric analysis and optimization in the Results section below.

The electrical analogy assumes linearized first-order governing equation between the volumetric flowrate of the working fluid, U_i through a given component i and the pressure difference, ΔP_i across that component:

$$\begin{array}{*{20}l} \Delta P_{i} = U_{i}\,Z_{i}\, {,} \end{array} $$

(2)

with Z_i the impedance of component i. The impedance for each of the main components of the engine is expressed in terms of resistance, inductance and capacitance in the Laplace domain, using the principle of superposition to describe multiple effects occurring within the same component:

$$\begin{array}{*{20}l} & \text{Load \textcircled{2}} \quad & Z_{\mathrm{l}} =\ & R_{\mathrm{l}} + s\,L_{\mathrm{l}}\\ & \text{Feedback connection \textcircled{6}:} \quad & Z_{\mathrm{f}} =\ & R_{\mathrm{f}} + s\,L_{\mathrm{f}}\\ & \text{Power cylinder \textcircled{7}:} \quad & Z_{\mathrm{p}} =\ & s\,L_{\mathrm{p}} + \frac{1}{s\,C_{\mathrm{p}}}\\ & \text{Displacer cylinder \textcircled{8}:} \quad & Z_{\mathrm{d}} =\ & s\,L_{\mathrm{d}} + \frac{1}{s\,C_{\mathrm{d}}}\\ & \text{Adiabatic vapour region \textcircled{9}:} \quad & Z_{\text{ad}} = & \frac{1}{s\,C_{\text{ad}}} \end{array} $$

These R,L,C elements, as arranged in a complete NIFTE circuit, may be visualized in Fig. 2, for which further details will be provided below.

Fig. 2

Complete NIFTE circuit. R_i denotes a resistance, C_i a capacitance, L_i an inductance, P_i a pressure, and U_i a volumetric flowrate. The subscript _th denotes the thermal domain, _ad the adiabatic domain, _l the load, _p the power cylinder, _d the displacer (heat-exchanger) cylinder, and _f the combined effect of the feedback connection and/or valve. The components in red relate to the part of the circuit describing irreversible thermal losses, only relevant to the NTP-TL model

Thermal domain

The thermal domain, which comprises the hot and cold heat-exchangers ➃, ➄ within the displacer cylinder ➇, converts heat into pressurization and fluid flow. We adopt the usual assumptions that: (i) the spatial temperature profile along the vertical walls of the heat-exchangers is externally imposed as a boundary condition to the device and does not vary with time; and (ii) the convective heat transfer associated with phase change dominates over the forced convection taking place away from an active region near the vapor-liquid interface of the working fluid within which phase-change heat transfer occurs. Under these conditions, the heat transferred to the working fluid, $\dot {X}_{\text {th}}$ is given by:

$$\begin{array}{*{20}l} \dot{X}_{\text{th}} = T^{\circ}\dot{S}_{\text{th}} = h_{\text{th}}\,A_{\text{th}}\,\left[T_{\text{th}}(y)-T_{\text{ad}}\right], \end{array} $$

(3)

with $\dot {S}_{\text {th}}$, the associated rate of change of entropy; T^∘, the equilibrium (time-averaged) temperature; h_th, the (constant) phase-change convective heat transfer coefficient; A_th, the surface area over which the heat transfer takes place; T_ad, the temperature of the working fluid; and T_th, the heat-exchanger wall temperature experienced by the working fluid at the vapor-liquid interface ➂, y.

The connection between Eq. 3 and the fluid domain proceeds by: (i) relating the entropy rate of change, $\dot {S}_{\text {th}}$ to an equivalent volumetric flow-rate of vapor (due to phase change), U_th via

$$\begin{array}{*{20}l} \dot{S}_{\text{th}} = \rho_{\text{vap}}^{\circ}\,\Delta s_{\text{vap}}\,U_{\text{th}}, \end{array} $$

with ρvap∘, the equilibrium (time-averaged) density of the vapor, and Δs_vap, the specific entropy of vaporization; and (ii) relating temperatures to pressures as

$$\begin{array}{*{20}l} T_{\text{th}} =\ & \left(\frac{\mathrm{d}T}{\mathrm{d}P}\right)_{\text{sat}} P_{\text{th}}, & T_{\text{ad}} =\ & \left(\frac{\mathrm{d}T}{\mathrm{d}P}\right)_{\text{sat}} P_{\text{ad}}, \end{array} $$

with $\left (\frac {\mathrm {d}T}{\mathrm {d}P}\right)_{\text {sat}}$, the rate of change in saturation temperature of the working fluid with respect to the saturation pressure. Using the same electrical analogy as in Eq. 2, we obtain:

$$\begin{array}{*{20}l} P_{\text{th}}-P_{\text{ad}} = f(P_{\mathrm{d}}) - P_{\text{ad}} = U_{\text{th}}\,Z_{\text{th}}, \end{array} $$

(4)

where the thermal impedance Z_th consists of the thermal resistance between the working fluid and the walls of the heat-exchangers,

$$\begin{array}{*{20}l} & \text{Displacer cylinder \textcircled{8}:} \quad & Z_{\text{th}} = R_{\text{th}} = \frac{\rho_{\text{vap}}^{\circ}\,\Delta s_{\text{vap}} T^{\circ}}{h_{\text{th}}\, A_{\text{th}} \left(\frac{\mathrm{d}T}{\mathrm{d}P}\right)_{\text{sat}}} \end{array} $$

The heat-exchanger wall temperature T_th, and the pressure P_th likewise, vary with the vertical position y of the vapor-liquid interface ➂, and can be related to the pressure variable P_d=ρ_wf g y. In particular, the function f in Eq. 4 may taken as either linear [6, 7] or nonlinear [25] depending on modeling assumptions.

In addition to the thermal resistance R_th, two important thermal loss mechanisms in the displacer cylinder ➇ are the shuttle loss and the entrance and retrograde condensation loss [6]. The former is associated with irreversible alternating conductive heat transfer between the working fluid and the cylinder walls, which may be modeled by a heat diffusion equation with either an isothermal or adiabatic boundary condition. The latter is caused by the alternating parasitic condensation and evaporation of the working fluid via convective heat transfer. Following the derivations in Refs. [8, 26], a total impedance representing these thermal losses can be solved in the Laplace domain, to arrive at the expressions:

$$\begin{array}{*{20}l} & \text{Isothermal boundary condition:}\\ & Z_{\text{tl}} = \frac{A_{\text{tl}}}{\sqrt{s}} \tanh\left(B_{\text{tl}}\sqrt{s}\right) + R_{\text{tl},0} \end{array} $$

(5a)

$$\begin{array}{*{20}l} & \text{Adiabatic boundary condition:}\\ & Z_{\text{tl}} = \frac{A_{\text{tl}}}{\sqrt{s}} \coth\left(B_{\text{tl}}\sqrt{s}\right) + R_{\text{tl},0} \end{array} $$

(5b)

where the thermal loss parameters A_tl,B_tl and R_tl,0 are given by

$${\begin{aligned} A_{\text{tl}} =\ & \frac{\rho_{\text{vap}}^{\circ} T^{\circ} \Delta s_{\text{vap}} \sqrt{\alpha_{\text{gl}}}}{\Sigma_{\text{gl}} \left(\frac{\mathrm{d}T}{\mathrm{d}P}\right)_{\text{sat}} \kappa_{\text{gl}}}, & B_{\text{tl}} =\ & \frac{l_{\text{gl}}}{\sqrt{\alpha_{\text{gl}}}}, & R_{\text{tl},0} =\ & \frac{\rho_{\text{vap}}^{\circ} T^{\circ} \Delta s_{\text{vap}}}{h_{\text{th}} \Sigma_{\text{gl}} \left(\frac{\mathrm{d}T}{\mathrm{d}P}\right)_{\text{sat}}} \end{aligned}} $$

with α_gl,κ_gl,l_gl and Σ_gl, the thermal diffusivity, thermal conductivity, length and surface area of the glass component, respectively. We reiterate that the model assumes the physical properties to be temperature independent.

In order to facilitate expressing the dynamic system in the time domain, we consider an RC-circuit approximation for the nonlinear expression of Z_tl herein, in the form:

$$\begin{array}{*{20}l} Z_{\text{tl}} \approx R_{\text{tl},0} + \sum\limits_{i=0}^{n} \left[ \frac{1}{R_{\text{tl},i}} + C_{\text{tl},i} s \right]^{-1}, \end{array} $$

(6)

where n denotes the number of parallel R-C components in series (see RC circuit depicted in red on Fig. 2). The resistance and capacitance values R_tl,i and C_tl,i may be estimated to match Z_tl over a suitable frequency range.

Complete NIFTE models

A graphical depiction of the R,L,C elements within each component of both the fluid and thermal domains is shown in Fig. 2. The original NTP model [25] comprises the components shown in black only. Addition of the components in red describing irreversible thermal losses, gives rise to the extended NTP-TL model.

A mathematical formulation of the NTP model (Eqs. 7a–7e) from its electrical circuit analogue is obtained by applying Kirchhoff’s voltage and current laws. A minimum of five independent thermodynamic variables is required to fully describe the system, here taken as: the pressure in the adiabatic vapor region, $\hat {P}_{\text {ad}}$; the hydrostatic pressures in the displacer cylinder, $\hat {P}_{\mathrm {d}}$ and in the power cylinders, $\hat {P}_{\mathrm {p}}$; and the volumetric flowrates in the feedback tube, $\hat {U}_{\mathrm {f}}$ and in the power cylinder $\hat {U}_{\mathrm {p}}$. The resulting system of ordinary differential equations (ODEs) is furthermore non-dimensionalized (denoted with $\hat {\,}$), by using P^∘,U^∘ and τ as reference values for the pressures, flowrates and time, respectively.

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{P}_{\text{ad}}}{\mathrm{d} \hat{t}} = & \frac{\tau\left[f\left(\hat{P}_{\mathrm{d}}\right)-\hat{P}_{\text{ad}}\right]}{R_{\text{th}}C_{\text{ad}}} + \frac{\tau U^{\circ}\left(\hat{U}_{\mathrm{f}}+\hat{U}_{\mathrm{p}}\right) }{P^{\circ} C_{\text{ad}}} \end{array} $$

(7a)

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{P}_{\mathrm{d}}}{\mathrm{d} \hat{t}} = & \frac{\tau U^{\circ} \hat{U}_{\mathrm{f}}}{P^{\circ} C_{\mathrm{d}}} \end{array} $$

(7b)

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{P}_{\mathrm{p}}}{\mathrm{d} \hat{t}} = & \frac{\tau U^{\circ} \hat{U}_{\mathrm{p}}}{P^{\circ} C_{\mathrm{p}}} \end{array} $$

(7c)

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{U}_{\mathrm{f}}}{\mathrm{d} \hat{t}} = & \frac{(\tau P^{\circ}/U^{\circ})\left[ L_{\mathrm{l}}\hat{P}_{\mathrm{p}} - L_{\mathrm{p}}\hat{P}_{\text{ad}} - (L_{\mathrm{p}}+L_{\mathrm{l}})\hat{P}_{\mathrm{d}} \right]}{L_{\mathrm{l}}(L_{\mathrm{d}}+L_{\mathrm{f}}) + L_{\mathrm{p}}(L_{\mathrm{d}}+L_{\mathrm{f}}+L_{\mathrm{l}}) } \end{array} $$

$$\begin{array}{*{20}l} & - \frac{\tau[ L_{\mathrm{l}} R_{\mathrm{f}} + L_{\mathrm{p}}(R_{\mathrm{f}} + R_{\mathrm{l}}) ]\hat{U}_{\mathrm{f}} + \tau L_{\mathrm{p}} R_{\mathrm{l}} \hat{U}_{\mathrm{p}}} {L_{\mathrm{l}}(L_{\mathrm{d}}+L_{\mathrm{f}}) + L_{\mathrm{p}}(L_{\mathrm{d}}+L_{\mathrm{f}}+L_{\mathrm{l}})} \\ \frac{\mathrm{d} \hat{U}_{\mathrm{p}}}{\mathrm{d} \hat{t}} = & \frac{(\tau P^{\circ}/U^{\circ})\left[ L_{\mathrm{l}}\hat{P}_{\mathrm{d}} \,-\, (L_{\mathrm{d}}+L_{\mathrm{f}})\hat{P}_{\text{ad}} \,-\, (L_{\mathrm{d}}+L_{\mathrm{f}}+L_{\mathrm{l}})\hat{P}_{\mathrm{p}} \right]}{L_{\mathrm{l}}(L_{\mathrm{d}}+L_{\mathrm{f}}) + L_{\mathrm{p}}(L_{\mathrm{d}}+L_{\mathrm{f}}+L_{\mathrm{l}})} \end{array} $$

(7d)

$$\begin{array}{*{20}l} & + \frac{\tau[ L_{\mathrm{l}} R_{\mathrm{f}} - (L_{\mathrm{d}}+L_{\mathrm{f}})R_{\mathrm{l}} ]\hat{U}_{\mathrm{f}} - \tau (L_{\mathrm{d}}+L_{\mathrm{f}})R_{\mathrm{l}} \hat{U}_{\mathrm{p}}} {L_{\mathrm{l}}(L_{\mathrm{d}}+L_{\mathrm{f}}) + L_{\mathrm{p}}(L_{\mathrm{d}}+L_{\mathrm{f}}+L_{\mathrm{l}})} \end{array} $$

(7e)

The derivation of the NTP-TL model (Eqs. 7a’–7f) proceeds similarly to that of the NTP model above, yet with n extra thermodynamic variables $\hat {P}_{\text {tl},i}$ representing the pressures associated with each thermal loss capacitance C_tl,i. Eq. 7a’ describing the dynamics of pressure $\hat {P}_{\text {ad}}$ also comprises an extra term compared with Eq. 7a, to account for the thermal losses.

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{P}_{\text{ad}}}{\mathrm{d} \hat{t}} = & \frac{\tau\left[f(\hat{P}_{\mathrm{d}})-\hat{P}_{\text{ad}}\right]}{R_{\text{th}}C_{\text{ad}}} + \frac{\tau U^{\circ}\left(\hat{U}_{\mathrm{f}}+\hat{U}_{\mathrm{p}}\right) }{P^{\circ} C_{\text{ad}}} - \frac{\tau\left(\hat{P}_{\text{ad}} - \sum\nolimits_{i=1}^{n} \hat{P}_{\text{tl},i}\right)}{R_{\text{tl},0}C_{\text{ad}}} \end{array} $$

(7a’)

$$\begin{array}{*{20}l} \frac{\mathrm{d} \hat{P}_{\text{tl},i}}{\mathrm{d} \hat{t}} =\ & \frac{\tau \left(\hat{P}_{\text{ad}} - \sum\nolimits_{i=1}^{n} \hat{P}_{\text{tl},i}\right)}{R_{\text{tl},0} C_{\text{tl},i}} - \frac{\tau \hat{P}_{\text{tl},i}}{R_{\text{tl},i} C_{\text{tl},i} }, \quad i=1\ldots n \end{array} $$

(7f)

Observe that both NTP and NTP-TL models consist of linear ODEs in the state variables, with the exception of the term $f(\hat {P}_{\mathrm {d}})$ in Eqs. 7a and 7a’ describing the transfer function between $\hat {P}_{\mathrm {d}}$ and $\hat {P}_{\text {th}}$. Early models of the NIFTE [4–7, 19] have considered linear (static) transfer functions of the form $\hat {P}_{\text {th}}=K_{\text {th}}\,\hat {P}_{\mathrm {d}}$, where K_th is a proportional gain. An alternative approach that also preserves linearity uses the (dynamic) transfer function $\hat {P}_{\text {th}}\,s=K_{\text {th}}\,\hat {P}_{\mathrm {d}}$ [8, 27]. But because they are fully linear, none of these models can capture the robust periodic oscillations corresponding to a (stable) limit cycle nonetheless.

To account for the fact that the temperature on the heat-exchanger walls T_th cannot increase or decrease indefinitely with the liquid level ➂ in the heat-exchanger block, one can introduce a nonlinearity in the form of a saturation function as the liquid level in the displacer cylinder moves away from the equilibrium position (which lies halfway between HHX ➃ and CHX ➄). Markides and co-workers [25] proposed the saturation function

$$\begin{array}{*{20}l} f(\hat{P}_{\mathrm{d}}) = \hat{P}^{\text{max}}_{\text{th}} \tanh\left(K_{\text{th}}\hat{P}_{\mathrm{d}}\right), \end{array} $$

(8)

where the parameter $\hat {P}^{\text {max}}_{\text {th}}$ sets the maximum (saturation) temperature amplitude on the heat-exchanger walls; and the nonlinear gain K_th determines the maximum rate of heat transfer and is a function of the physical distance between the hot and cold heat elements. More specifically, K_th is inversely proportional to the distance between the HHX and the CHX; and directly proportional to the temperature difference between the heat source and sink. The resulting NIFTE models are comprised of nonlinear ODEs, which have been shown to feature limit cycles [25].

Key performance indicators

An important indicator of the performance of a NIFTE is the average useful power output from the load, $\dot {W}_{\mathrm {l}}$ given by

$$\begin{array}{*{20}l} \dot{W}_{\mathrm{l}} = \overline{P_{\mathrm{l}}(t)U_{\mathrm{l}}(t)} \end{array} $$

(9)

where P_l and U_l denote the pressure and volumetric flowrate in the load, respectively; and the average value is calculated over a steady-state cycle of the engine.

A second key performance criterion is the exergetic (or, second law) efficiency of the engine, η_ex defined as

$$\begin{array}{*{20}l} \eta_{\text{ex}} = \frac{\dot{W}_{\mathrm{l}}}{\dot{X}_{\text{th}}} = \frac{\overline{ P_{\mathrm{l}}(t)U_{\mathrm{l}}(t)}} {\overline{P_{\text{th}}(t)U_{\text{th}}(t)}} \end{array} $$

(10)

where P_th and U_th are the input thermal pressure and volumetric flowrate associated with the phase change, respectively; and $\dot {X}_{\text {th}}$ the exergy power input into the device.

For given temperatures of the heat source and sink, recall that the Carnot efficiency, η_ca is constant, and therefore η_ex is proportional to the thermal efficiency of the engine, η_th=η_ex η_ca. The extra pressure and volumetric flowrate states in Eqs. 9–10 are readily determined from the application of Kirchoff’s circuit laws in the electrical circuit analogue (Fig. 2). Further information on the definition and interpretation of these performance indicators may be found in Reference [7].

Methods

Cyclic steady-state detection and analysis

The search for a CSS of the NIFTE can be automated by solving the following constrained dynamic optimization problem,

$$\begin{array}{*{20}l} \min_{{\mathbf{x}_{0}, T}}\ & T \\ \textrm{s.t.}\ \ & \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}) \\ & \mathbf{x}(0) = \mathbf{x}(T) = \mathbf{x}_{0} \end{array} $$

(11)

where T denotes a cycle period; the ODEs right-hand side, f corresponds to either the NTP model (Eqs. 7a–7e) with the state vector x comprised of $\left (\hat {P}_{\text {ad}}, \hat {P}_{\mathrm {d}}, \hat {P}_{\mathrm {p}}, \hat {U}_{\mathrm {f}}, \hat {U}_{\mathrm {p}}\right)$, or the NTP-TL model (Eqs. 7a’–7f) with x comprised of $\left (\hat {P}_{\text {ad}}, \hat {P}_{\mathrm {d}}, \hat {P}_{\mathrm {p}}, \hat {P}_{\text {tl}}, \hat {U}_{\mathrm {f}}, \hat {U}_{\mathrm {p}}\right)$; and p is the vector of model parameters (cf. Table 2).

Table 2 Nominal parameter values in NTP and NTP-TL models of the NIFTE

Any solution $\overline {\mathbf {x}}(\cdot,\mathbf {p})$ to the problem (11) with T>0 corresponds to a CSS of the NIFTE. However, this optimization problem is over-parameterized, in the sense that there exists infinitely many initial conditions x₀ describing a given CSS with T>0. This issue is easily circumvented – and without loss of generality – by forcing one of the initial states to zero; e.g., $\hat {P}_{\text {ad}}(0)=0$.

The optimization problem (11) is nonconvex and may therefore present multiple (local) optima. These optima may correspond to different CSS or equilibrium points of the engine, but they could also be “parasitic” local optima corresponding to multiple periods, e.g. 2T,3T,…, of the same CSS. Herein, we use a simple single-shooting approach as implemented in our in-house optimization toolkit CRONOS¹ to compute such local optima, and we apply a multi-start search heuristic based on Sobol’ sampling to detect the presence of multiple local optima. CRONOS provides an interface to the nonlinear programming (NLP) solver IPOPT [28] coupled with the ODE integrator CVODES in SUNDIALS [29] for computing the objective/constraint functions and their gradients. A local optimization typically converges within a few seconds of CPU time.

The stability of the limit cycle corresponding to a solution $\overline {\mathbf {x}}(\cdot,\mathbf {p})$ can be characterized by analyzing the so-called monodromy matrix Φ_p(T,0), which is the fundamental matrix solution to the variational equation [30, 31]

$$\begin{array}{*{20}l} & \forall \tau, t\in[0,T],\ \ \frac{\partial \Phi_{\mathbf{p}}(t,\tau)} {\partial t} = \frac{\partial \mathbf{f}(t, \overline{\mathbf{x}}(t,\mathbf{p}), \mathbf{p}) }{\partial \mathbf{x}} \Phi_{\mathbf{p}}(t, \tau) \end{array} $$

(12a)

$$\begin{array}{*{20}l} & \text{with}\ \ \Phi_{\mathbf{p}}(\tau,\tau) = I \end{array} $$

(12b)

In practice, the monodromy matrix can be conveniently computed as $\Phi _{\mathbf {p}}(T,0) = \frac {\partial \overline {\mathbf {x}}(T,\mathbf {p})}{\partial \mathbf {x}_{0}}$, e.g. using the sensitivity analysis capabilities of the integrator CVODES.

One (at least) of the eigenvalues of the monodromy matrix shall always be equal to 1. The limit cycle $\overline {\mathbf {x}}(\cdot,\mathbf {p})$ is asymptotically stable (attractor) if all of the remaining eigenvalues lie within the open unit disk,

$$ \max \limits_{i} | \lambda_{i} (\Phi_{\mathbf{p}}(T,0)) | < 1 $$

(13)

If any of these eigenvalues has a modulus greater than 1, the limit cycle is otherwise unstable (repeller).

Bifurcation occurs when a small change made to the values of a system’s parameters p causes a sudden qualitative change in its asymptotic behavior [31]. Of particular interest in the NIFTE are so-called Hopf bifurcations, whereat the equilibrium position x=0 (origin) loses stability and a periodic solution arises. Such bifurcations can be detected upon repeating the solution of the optimization problem (11) at a large number of parameter values. One can then use dedicated codes to confirm the presence of a bifurcation and refine its location, such as MatCont [32, 33] which implements a numerical continuation approach in MATLAB. The presence of multiple Hopf bifurcations for a particular parameter can be an indication for the co-existence of multiple limit cycles for certain engine configurations.

Multi-objective optimization

Finding the design parameters that maximize both the exergetic efficiency and the average power output of a NIFTE leads to a multi-objective optimization (MOO) problem in the form of

$$\begin{array}{*{20}l} \min_{\mathbf{x}_{0}, T, \mathbf{p}}\ & \{\dot{W}_{\mathrm{l}}(\mathbf{x}(\cdot)), \eta_{\text{ex}}(\mathbf{x}(\cdot)) \} \\ \textrm{s.t.}\ \ & \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}) \\ & \mathbf{x}(0) = \mathbf{x}(T) = \mathbf{x}_{0} \\ & \mathbf{p}^{\mathrm{L}} \leq \mathbf{p} \leq \mathbf{p}^{U} \end{array} $$

(14)

with $\dot {W}_{\mathrm {l}}$ and η_ex as defined in Eqs. 9–10; and using the same notation as in the problem (11). Notice that the initial conditions x₀ and cycle time T are decision variables alongside the design variables p in order to enforce optimal performance under CSS operation. The solution set to the problem (14) gives rise to a so-called Pareto optimal frontier [34], whereby each point is non-dominated in the sense that no better feasible solution exist in terms of the two objectives simultaneously: betterment of the exergetic efficiency compared to a point on the Pareto frontier results in worsening of the average power output, and vice versa.

Herein, we approximate the Pareto frontier of the problem (14) via scalarization, and we proceed by solving the following single-objective optimization problems for a collection of scalars α∈[0,1]:

$$\begin{array}{*{20}l} \min_{\mathbf{x}_{0}, T, \mathbf{p}}\ & (1-\alpha) \dot{W}_{\mathrm{l}}(\mathbf{x}(\cdot)) + \alpha \eta_{\text{ex}}(\mathbf{x}(\cdot)) \\ \textrm{s.t.}\ \ & \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}) \\ & \mathbf{x}(0) = \mathbf{x}(T) = \mathbf{x}_{0} \\ & \mathbf{p}^{\mathrm{L}} \leq \mathbf{p} \leq \mathbf{p}^{\mathrm{U}} \end{array} $$

(15)

Clearly, the problems with α=0 and α=1 correspond to the two extremes of maximizing the average power output and the exergetic efficiency, respectively. It is noteworthy that such a simple scalarization approach fails to describe nonconvex portions of a Pareto frontier in general [34], but our results later on will confirm that this situation does not occur for the Pareto frontiers of interest herein.

In order to tackle the single-objective optimization problem (15), we furthermore apply a two-level decomposition: for a given design vector p, the inner problem determines the values of x₀ and T that yield the CSS with best performance, e.g. in case multiple CSS coexist, or detects that no (nontrivial) CSS is possible; whereas the outer optimization problem operates on the design vector p only, and solves the inner optimization problem at each iteration in order to calculate the performance under CSS operation for a particular instance of p. The inner optimization problems are solved in the same way as the CSS detection problem (11) earlier, i.e. using a combination of local search and multi-start search heuristic. As for the outer problem, we implement a simple Nelder-Mead simplex algorithm [35], a derivative-free approach with slow convergence but tolerant to numerical noise, with a termination tolerance on the minimum simplex diameter (set to 10⁻³ in log-scale here). In order to rule out convergence to local (suboptimal) solutions and initialize the Nelder-Mead simplex algorithm, the search space is probed using low-discrepancy sampling based on Sobol’ sequences [36] (with 10,000 points here), and the NIFTE performance under CSS operation is computed at each sampling point. Algorithm 1 provides a summary of the multi-objective optimization procedure. The parameter range is set to ± 3 orders of magnitude around their nominal values (see below), and the set of weights considered in the scalarized cost of the problem (15) are $\left \{0,\frac {1}{2000},\frac {1}{1000},\frac {1}{500},\frac {1}{200},\frac {1}{100},\frac {1}{50},\frac {1}{20},\frac {1}{10},\frac {1}{5},1\right \}$. Each run of the Nelder-Mead simplex algorithm in order to determine one point on the Pareto frontier takes up to a few hundred iterations to converge.

Results and discussion

Case study description

The assessment presented in this section is based on a working prototype of the NIFTE. A temperature difference of 80 ^∘C is applied across the heat-exchanger, with the cold source at ambient temperature and a corresponding Carnot efficiency of about 20%. The pumped and working fluids in this prototype are water and n-pentane, respectively. The load corresponds to a delivery pipe of 1 m in length and with a diameter of 0.01 m that would in practice be connecting a supply tank of water to the destination tank. Further details concerning this prototype may be found in Refs. [6, 7, 19].

Table 2 reports the nominal values for the electrical analogues R_i,L_i and C_i in the NTP model (Eqs. 7a–7e) as well as the parameters $K_{\text {th}}, \hat {P}_{\text {th}}^{\text {max}}$ in the heat transfer function (Eq. 8) and the reference values P^∘,U^∘,τ. These were calculated based on the prototype’s physical properties, and it is noteworthy that the NTP model was tested against experimental data with these parameter values (referred to as Set-I in Ref. [25]). In spite of its simplicity compared to complexity of the actual device, the NTP model is capable of early-stage engineering predictions of the NIFTE’s operation and performance.

The same parameter values are used in the NTP-TL model (Eqs. 7a’–7f). Values for the extra parameters n, R_tl,i,C_tl,i in the approximate thermal-loss impedance (Eq. 6) were determined by fitting the thermal-loss impedance with isothermal boundary condition (Eq. 5a), using standard MATLAB tools. In order to obtain accurate gain and phase predictions in the frequency range of interest, 0.1– 100 rad s⁻¹ for the NIFTE prototype at hand, a minimum of n=4 RC-circuits in series is necessary here (Fig. 3). The fitted values of R_tl,i and C_tl,i for i=1…4 are reported in Table 2. Notice that the values of the parameters A_tl,B_tl,R_tl,0 in Eq. 5 used for the calibration were also calculated based on the prototype’s physical properties, reported in Table 2 for completeness. Interestingly, the use of RC-circuit in series to approximate the thermal-loss impedance with isothermal boundary condition also yields an excellent agreement for its counterpart with adiabatic boundary condition (Eq. 5b) in the frequency range 0.2– 100 rad s⁻¹. Such close agreement suggests that both types of boundary conditions in Eq. 5 to express the thermal-loss impedance should lead to minimal differences in quantitative predictions by the NTP-TL model. Both types of boundary conditions are indeed relevant as they are representative of different applications of the NIFTE technology.

Fig. 3

Bode plot of thermal loss impedance. The approximate thermal-loss impedance (Eq. 6) with n=4 and values of R_tl,i,C_tl,i as in Table 2 is plotted in black solid line. The thermal-loss impedances with isothermal boundary condition (Eq. 5a) and adiabatic boundary condition (Eq. 5b) are plotted in blue and red dashed lines, respectively

Cyclic-steady-state multiplicity and stability

This subsection investigates the asymptotic behavior of a NIFTE based on the NTP and NTP-TL models. It has already been established [25], based on the NTP model, that a given engine may or may not exhibit a limit cycle. A key finding herein is that a given engine could in fact exhibit multiple CSS, which could be either stable or unstable, an observation that had never been reported before.

To illustrate this behavior, we first consider a NIFTE with all its parameters at their nominal values (Table 2), except the feedback inductance L_f=4.74×10⁵ kg m⁻⁴ s⁻¹ (10-fold smaller than the nominal value). When using the NTP model and solving the optimization problem (11) as described earlier, three possible asymptotic behaviors are identified for this NIFTE configuration. The first one is a trivial equilibrium point at the origin, which is unstable. The other two are limit cycles with periods T=1.30 s (CSS1) and T=0.88 s (CSS2). The spectra of the monodromy matrices (Eq. 12) at CSS1 and CSS2 are reported in Table 3. Since all the eigenvalues of the monodromy matrix at CSS1 lie within the unit disk (Eq. 13), this CSS describes a stable limit cycle. By contrast, the monodromy matrix at CSS2 has two of its eigenvalues outside the unit disk (shown in red), and so this is an unstable limit cycle.

Table 3 Asymptotic stability and performance of a NIFTE with nominal parameter values (Table 2), except L_f=4.74×10⁵ kg m⁻⁴ s⁻¹

The pressure-volume phase diagrams in Fig. 4 are obtained by simulating the NIFTE over an extended time horizon, from different initial states. They represent the pressure P_l against the volume V_l of the working fluid, given by $\dot {V}_{\mathrm {l}}=U_{\mathrm {l}}$, and so the perimeter of the cycle limit is proportional to the work output, W_l over one complete cycle:

$$\begin{array}{*{20}l} W_{\mathrm{l}} = \int_{0}^{T} P_{\mathrm{l}}(t)\,U_{\mathrm{l}}(t)dt = \oint P_{\mathrm{l}}\,\mathrm{d}V_{\mathrm{l}} \end{array} $$

Fig. 4

Pressure-volume phase diagrams of a NIFTE simulated over an extended time horizon. The NIFTE is simulated using either the NTP or NTP-TL model with nominal parameters values (Table 2), except L_f=4.74×10⁵ kg m⁻⁴ s⁻¹. a: Trajectory of the NTP model repelled from the origin and attracted to CSS1. b: Trajectory of the NTP model repelled by CSS2 and attracted to CSS1. c: Trajectory of the NTP-TL model repelled from the origin and attracted to CSS1’

It is clear that CSS1 acts as an attractor for the trajectories, and the origin as a repeller (Fig. 4a). Furthermore, trajectories initialized near CSS2 eventually escape away from that limit cycle and are again attracted by CSS1, thereby showing two clusters of trajectories (Fig. 4b). All of these observations are in agreement with the numerical stability results in Table 3.

The comparison of performance indicators in the bottom half of Table 3 reveals that a 70%-greater average useful power $\dot {W}_{\mathrm {l}}$ (Eq. 9) is generated by operating the NIFTE at the (unstable) CSS2 compared with the (stable) CSS1. Because the engine is attracted to CSS1 when initialized near the origin (Fig. 4a), reaching (and then remaining around) an alternative, potentially better limit cycle such as CSS2 would require a tailored control strategy. This could be done, for example, by manipulating the valve in the feedback connection ➅ (Fig. 1), a dynamic control that has never been applied in practice however. Notice also that CSS2 has a 30%-shorter period T and an 8-fold greater exergetic efficiency η_ex (Eq. 10) than CSS1. The latter is not only the result of the greater useful power output at CSS2, but also a 5-fold lower exergy power input $\dot {X}_{\text {th}}$ (Eq. 10) into the engine from the heat source. Coincidentally, the combination of a shorter period and a larger power output at CSS2 results in about the same work output during a complete cycle of CSS1 (W_l=1.03 J) or CSS2 (W_l=1.07 J). This can seen on Fig. 4b, where the limit cycles for CSS1 and CSS2 appear to have identical perimeters.

Unlike the NTP model, a simulation of the same NIFTE with the NTP-TL model in problem (11) predicts a unique, stable limit cycle (CSS1’). This CSS has a slightly longer period than CSS1, while both the exergetic efficiency and power output are reduced by about 3-fold (Table 3). This dramatic reduction is of course attributed to the extra thermal losses in NTP-TL, confirming that their inclusion is paramount to the accuracy of the model predictions. The perimeter of the limit cycle CSS1’ on the pressure-volume phase diagrams in Fig. 4c is accordingly much shorter in comparison with that of CSS1 on Fig. 4a.

In the following subsections, a parametric analysis of the NIFTE is carried out by varying (combinations of) the design parameters in the feedback, displacer and power cylinders and heat-exchanger, and analyzing the performance of the corresponding CSS using both the NTP and NTP-TL models. A multi-objective optimization is also conducted with the NTP-TL model in order to investigate the optimal trade-offs between exergetic efficiency and power input when varying three or more parameters.

Parametric analysis of (R _f,L _f)

Previous studies [6, 7] have shown that the performance of NIFTE is highly sensitive to the feedback resistance R_f and the feedback inductance L_f. These parameters can be varied in a physical device by modifying the geometry of the feedback tube (e.g., l_f,d_f or A_f; see Eq. 1) or the position of a feedback valve thereof. This subsection presents a joint parametric analysis of the feedback parameters R_f and L_f using both models, while all other parameters are held constant at their nominal values (Table 2). We use the superscript ⁺ to denote the ratio between a parameter and its nominal value.

NTP model

Both the exergetic performance η_ex and useful power output $\dot {W}_{\mathrm {l}}$ predicted by the NTP model under CSS operation are shown in Fig. 5, for a wide parameter range spanning several orders of magnitude around the nominal parameter values. By repeating the solution of the optimization problem (11) over this range, we could identify up to two limit cycles, denoted by CSS1 and CSS2 in agreement with previous notation. The stability of each of these CSS is further determined by the eigenvalues of their monodromy matrices (Eq. 12), and we report the separation between stable and unstable CSS on each plot in Fig. 5. Notice that multiple CSS only co-exist in the subdomain $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right)\in [0.001,10]\times [0.044,0.145]$, with either one stable and one unstable CSS, or even two stable CSS.

Fig. 5

Predicted effect of R_f and L_f on NIFTE performance using the NTP model. Two distinct cyclic steady states are detected, denoted by CSS1 and CSS2, which are overlapping for certain (R_f,L_f) combinations. The effect on both the exegetic efficiency η_ex (left plots) and the average power output $\dot {W}_{\mathrm {l}}$ (right plots) is shown for each cyclic steady-state. $R_{\mathrm {f}}^{+}$ and $L_{\mathrm {f}}^{+}$ refer to R_f and L_f normalized by their respective nominal values in Table 2. The remaining parameters are set to their nominal values. The separatrix between stable and unstable limit cycles are shown in white, with arrows pointing in the direction of stable limit cycles

As far as the effect of R_f is concerned, a reduction of the feedback resistance with all other parameters fixed leads to a larger useful power output and a larger exergetic efficiency, for both CSS1 (Fig. 5a) and CSS2 (Fig. 5b). This behavior is consistent with previous parametric analyses conducted with the NTP model [25]. The parameter L_f is also found to have a high impact on the exergetic performance η_ex and useful power output $\dot {W}_{\mathrm {l}}$, but analyzing these variations is more intricate. In the case of CSS1, an overall trend is for the useful power output to increase as L_f is lowered, whereas the exergetic efficiency is maximal around the nominal L_f value $\left (L_{\mathrm {f}}^{+}=1\right)$. These different trends, alongside the wavy pattern on Fig. 5a, is attributed to a shift in phase between oscillations in the load and displacer cylinders, as a result of the variation in gas-liquid interface and hydraulic pressure due to the phase change. Notice also that nearly all of CSS1 points are stable, apart from those corresponding to the smallest L_f values that turn into unstable before vanishing. In the case of CSS2 by contrast, a trend is for the exergetic efficiency to decrease as L_f is lowered, while the useful power output passes through a maximum that is shifting to lower L_f values as R_f itself is lowered.

The results in Fig. 5 thus reveal a trade-off between η_ex and $\dot {W}_{\mathrm {l}}$ in terms of maximizing NIFTE performance. The maximal efficiency η_ex is predicted to be as high as 77% at $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right) = (0.001,0.145)$, yet with a very small power output of 0.1 W; whereas the maximum power output is 5.4 W at $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right) = (0.001,0.001)$, but with a lower exergetic efficiency of about 36%. These values should be considered relative to the nominal NIFTE configuration $\left (R_{\mathrm {f}}^{+}=L_{\mathrm {f}}^{+}=1\right)$, which has a predicted exergetic efficiency and average power output of about 2.6% and 0.86 W, respectively. Whether the primary target is to maximize the efficiency or to produce more power should normally drive the design of a NIFTE in practice. Here, a good compromise could aim for a design near the separatrix between stable and unstable CSS2 in Fig. 5b; for instance at $L_{\mathrm {f}}^{+} = 0.03$ and $R_{\mathrm {f}}^{+}<0.01$, where the predicted power output and exergetic efficiency are greater than 5.2 W and 45%, respectively. Because these CSS are at the verge of instability, a control strategy [37] might become necessary in order to stabilize the NIFTE nonetheless.

Although the NTP model is capable of accurate predictions of several key operational characteristics of the NIFTE, such as its oscillation frequency, we shall see in the remainder of the paper that such improvements in efficiency and power output are unfortunately an artifact of this model as it fails to describe irreversible thermal losses.

NTP-TL model

We conduct the same parametric analysis of NIFTE, now accounting for thermal losses with the NTP-TL model. The predicted exergetic efficiency η_ex and average power output $\dot {W}_{\mathrm {l}}$ under CSS operation are presented in Fig. 6, for a wide parameter range spanning several orders of magnitude around the nominal parameter values. A first notable difference with the NTP-TL model is that a unique limit cycle is detected via the repeated solution of the optimization problem (11) for $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right) \in [\!0.001,10]^{2}$, which we refer to as CSS1’ in agreement with Table 3. This CSS1’ is furthermore stable over the explored parameter range.

Fig. 6

Predicted effect of R_f and L_f on NIFTE performance using the NTP-TL model. A unique, stable cyclic steady state is detected here. The effect on both the exegetic efficiency η_ex (left plots) and the average power output $\dot {W}_{\mathrm {l}}$ (right plots) is shown for each cyclic steady state. $R_{\mathrm {f}}^{+}$ and $L_{\mathrm {f}}^{+}$ refer to R_f and L_f normalized by their respective nominal values in Table 2. The remaining parameters are set to their nominal values. All of the limit cycles are stable on both plots

The variations in exergetic efficiency and power output of CSS1’ do not display the same wavy pattern that could be observed for CSS1 in Fig. 5a either, although CSS1’ presents several minima and maxima across the parameter range. The maximal efficiency η_ex is predicted to be around 1.1% at $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right)=(3.16,0.001)$, while a maximal average power output $\dot {W}_{\mathrm {l}}$ of 0.31 W is predicted at $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right)=(1.38,0.316)$. Thee values are now to be compared with the predicted exergetic efficiency of 0.69% and predicted average power output of 0.24 W for the nominal NIFTE configuration $\left (R_{\mathrm {f}}^{+}=L_{\mathrm {f}}^{+}=1\right)$. A general trend on Fig. 6 is for both the efficiency and power output to increase as L_f is reduced from its nominal value. In particular, these results suggest that enhanced performance of the NIFTE prototype by 30–50% could be achieved via a significant reduction in L_f, while maintaining R_f close to nominal.

These predictions by the NTP-TL model are a major departure from those of the NTP model. The comparison between Figs. 5a and 6a, where the feasible ranges of (R_f,L_f) largely overlap, shows that accounting for thermal losses may cause a 10- to 20-fold drop in the predicted exergetic efficiency and average power output, especially for the most promising NIFTE configurations determined by the NTP model. Despite the fact that the NTP and NTP-TL model predictions are close for certain operational characteristics of the NIFTE, such as its oscillation frequency, these extreme differences in efficiency and power output predictions shed doubt on whether the NTP model is capable of accurate, quantitative predictions of the NIFTE performance over wide parameter ranges. By contrast, the level of performance predicted by the NTP-TL model is in good agreement with direct experimental evidence [6, 7], where the exergetic efficiency of the NIFTE prototype was measured between 0.4–1.6%. It was also observed experimentally that the exergetic efficiency could improve by increasing the feedback resistance R_f, a behavior that is not described by the NTP model.

Parametric analysis of (C _p,L _p)

Apart from the feedback cylinder, the design of the power cylinder where the working fluid is in direct contact with the pumped fluid in the load (Fig. 1), is also known to impact the NIFTE performance greatly. This subsection carries out a joint parametric analysis of the capacitance C_p and inductance L_p based on the NTP-TL model only, with all the other parameters held constant at their nominal values (Table 2). Both C_p and L_p could be varied in a physical device by modifying either the length l_p or the cross-section area A_p of the power cylinder (Eq. 1).

The predicted exergetic efficiency η_ex and average power output $\dot {W}_{\mathrm {l}}$ under CSS operation are presented in Fig. 7, for $\left (L_{\mathrm {p}}^{+},C_{\mathrm {p}}^{+}\right) \in [0.01,100]^{2}$. By repeating the solution of the optimization problem (11) over this parameter range, we could identify three distinct limit cycles, denoted by CSS1”, CSS2” and CSS3”. Notice that CSS1” is the unique limit cycle in a wide domain around the nominal NIFTE design $\left (L_{\mathrm {p}}^{+} = C_{\mathrm {p}}^{+} = 1\right)$. The domains of existence of CSS2” and CSS3” overlap with CSS1” for values of $C_{\mathrm {p}}^{+}\leq 0.025$ and $L_{\mathrm {p}}^{+}\geq 17$, respectively, yet these domains are much narrower.

Fig. 7

Predicted effect of C_p and L_p on NIFTE performance using the NTP-TL model. Three distinct cyclic steady states are detected, denoted by CSS1”, CSS2”, and CSS3”, which are overlapping for certain (C_p,L_p) combinations. The effect on both the exegetic efficiency η_ex (left plots) and the average power output $\dot {W}_{\mathrm {l}}$ (right plots) is shown for each cyclic steady-state. $C_{\mathrm {p}}^{+}$ and $L_{\mathrm {p}}^{+}$ refer to C_p and L_p normalized by their respective nominal values in Table 2. The remaining parameters are set to their nominal values. The separatrix between stable and unstable limit cycles are shown in white, with arrows pointing in the direction of stable limit cycles. All of the limit cycles are stable on those plots that do not have a separatrix

By far the most promising limit cycle is CSS1” in Fig. 7a, both in terms of exergetic efficiency and power output. The entire range of CSS1” is furthermore stable. The maximal efficiency η_ex is predicted to be slightly below 1% at $\left (C_{\mathrm {p}}^{+},L_{\mathrm {p}}^{+}\right) = (24.5, 0.105)$, with efficiencies above 0.9% achieved in a wide subdomain of CSS1” with $C_{\mathrm {p}}^{+}\geq 3$ and $L_{\mathrm {p}}^{+}\leq 1$. A maximal average power output $\dot {W}_{\mathrm {l}}$ of 0.31 W is predicted at $\left (C_{\mathrm {p}}^{+},L_{\mathrm {p}}^{+}\right) = (4.90,0.01)$, with improved power outputs again observed in a wide band with $3\leq C_{\mathrm {p}}^{+}\leq 10$ and $L_{\mathrm {p}}^{+}$ around or below nominal. Similar to the sensitivity analysis for (R_f,L_f) earlier, these results suggest that the performance of the NIFTE prototype could be enhanced by 30–50% by increasing C_p and reducing L_p from their nominal values in the NIFTE prototype at hand.

Parametric analysis of (C _d,L _d)

A third sensitivity analysis is conducted using the NTP-TL model in terms of the capacitance C_d and inductance L_d of the displacer cylinder, with all the other parameters held constant at their nominal values (Table 2). Similar to the parameter pair (C_p,L_p) earlier, both C_d and L_d could be varied in a physical device by modifying either the length l_d or the cross-section area A_d of the displacer cylinder (Eq. 1).

The predicted exergetic efficiency η_ex and average power output $\dot {W}_{\mathrm {l}}$ under CSS operation are shown in Fig. 8, for $\left (L_{\mathrm {d}}^{+},C_{\mathrm {d}}^{+}\right) \in [0.01,100]^{2}$. By repeating the solution of the optimization problem (11) over this parameter range, two distinct limit cycles emerged, denoted by CSS1”’ and CSS2”’. Like previously with (C_p,L_p), we find that CSS1”’ is the unique limit cycle in a wide subdomain around the nominal NIFTE design $\left (L_{\mathrm {d}}^{+} = C_{\mathrm {d}}^{+} = 1\right)$. The domain of existence of CSS2”’ overlaps with that of CSS1”’ over a narrow band with $C_{\mathrm {d}}^{+} \in [0.046,0.063]$. A majority of the limit cycles are furthermore stable, under either CSS1”’ or CSS2”’.

Fig. 8

Predicted effect of C_d and L_d on NIFTE performance using the NTP-TL model. Two distinct cyclic steady states are detected, denoted by CSS1”’ and CSS2”’, which are overlapping for certain (C_d,L_d) combinations. The effect on both the exegetic efficiency η_ex (left plots) and the average power output $\dot {W}_{\mathrm {l}}$ (right plots) is shown for each cyclic steady-state. $C_{\mathrm {d}}^{+}$ and $L_{\mathrm {d}}^{+}$ refer to C_d and L_d normalized by their respective nominal values in Table 2. The remaining parameters are set to their nominal values. The separatrix between stable and unstable limit cycles are shown in white, with arrows pointing in the direction of stable limit cycles

The most promising limit cycle is CSS1”’ in Fig. 8a, with predicted exergetic efficiencies and power outputs three- to four-times higher than CSS2”’ in Fig. 8b. A maximal efficiency η_ex of about 0.93% is predicted at $\left (C_{\mathrm {d}}^{+}, L_{\mathrm {d}}^{+}\right) = (0.229, 0.010)$; and the maximum power output $\dot {W}_{\mathrm {l}}$ is predicted to be 0.32 W at $\left (C_{\mathrm {d}}^{+}, L_{\mathrm {d}}^{+}\right) = (0.055, 95.5)$. The gain in power output is about identical to the gains predicted by adjusting the feedback parameters (R_f,L_f) or the displacer parameters (C_p,L_p); whereas the gain in efficiency, in the order of 35% compared to the nominal NIFTE, is lower than the gains predicted by adjusting by (R_f,L_f) or (C_p,L_p). A general trend on Fig. 8a is for both the efficiency and power output to increase as C_d is reduced to within the range $C_{\mathrm {d}}^{+}\in [0.15,0.3]$, while the parameter L_d exhibits a relatively lower sensitivity in the NIFTE prototype.

Parametric analysis of K _th

To complete the assessment, we conduct a final parametric analysis of the nonlinear gain K_th using the NTP-TL model, with all the other parameters held constant at their nominal values (Table 2). The gain represents the maximum rate of heat transfer across the heat-exchanger, which may be varied in a physical device by adjusting either the distance between the HHX and CHX elements (Fig. 1), or the temperature difference thereof.

The top-left plot in Fig. 9 shows the predicted exergetic efficiency η_ex versus average power output $\dot {W}_{\mathrm {l}}$ under CSS operation, for $K_{\text {th}}^{+} \in [0.09, 10]$. The embedded table reports the values of $K_{\text {th}}^{+}$ at four selected points (red circles), in addition to the nominal point (green triangle). The general trend is for η_ex to increase and $\dot {W}_{\mathrm {l}}$ to decreases as K_th is reduced, i.e. as the HHX and CHX elements are drawn further apart or the temperature difference is reduced. The sensitivity of $K_{\text {th}}^{+}$ appears to be very low around the nominal value $\left (K_{\text {th}}^{+}=1\right)$, where increasing $K_{\text {th}}^{+}$ only leads to marginal variations in average power output and exegetic efficiency. The engine performance becomes notably more sensitive to values of $K_{\text {th}}^{+}<0.5$, but a CSS ceases to exist below $K_{\text {th}}^{+} \approx 0.09$ due to the HHX and CHX elements no longer allowing for a sufficient heat transfer rate. The exergetic efficiency there is predicted to increase by over 10% compared to the nominal efficiency, but the NIFTE may no longer be practical because of the dramatic reduction in power output that ensues.

Fig. 9

Predicted Pareto frontiers for various parameter subsets of the NTP-TL model. The Pareto frontiers are for the exegetic efficiency η_ex against the average power output $\dot {W}_{\mathrm {l}}$. The subplots correspond to varying the parameters K_th (top left), R_f,L_f (top right), R_f,L_f,C_p (bottom left), and R_f,L_f,C_p,C_d (bottom right), with the other parameters kept at their nominal values in Table 2. Optimized (normalized) parameter values corresponding to the labeled points on the Pareto frontiers are reported in the tables on each subplot. The reference NIFTE design is marked with a green triangle. The CSS performance at sampled design parameters using Sobol’ sequences is displayed with orange dots

A physical interpretation of this parametric analysis is that the performance of the NIFTE would be seldom affected if either the temperature difference or the distance between heat-exchangers were to vary by a factor of 2 or even larger. In particular, these results corroborate experimental evidence and suggest that a NIFTE might indeed remain feasible for temperature differences lower than 30 ^∘C, without reducing the power output or exergetic efficiency significantly.

Multi-objective optimization by varying multiple parameters

Results of the parametric analyses conducted in the previous subsections revealed that adjusting the pairs (R_f,L_f),(C_p,L_p) and (C_d,L_d) of design parameters in the different NIFTE cylinders could enhance the work output by up to 30% or the exergetic efficiency by over 50% in the NIFTE prototype. For instance, the top-right plot in Fig. 9 shows the Pareto frontier of maximal exergetic efficiency versus maximal power output for the pair (R_f,L_f). This Pareto frontier was constructed using the same parametric sensitivity data as in Fig. 6. Values of $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+}\right)$ corresponding to the labeled points A–E on the frontier are also reported in the embedded table for comparison. Moving along the Pareto frontier from A to E in this case entails increasing the value of $R_{\mathrm {f}}^{+}$ by 2-fold, while simultaneously decreasing the value of $L_{\mathrm {f}}^{+}$ toward zero. Such Pareto frontiers are helpful to visualize the trade-offs between objectives, identify an acceptable compromise based on engineering priorities, and assess the sensitivity of Pareto-optimal decisions toward the design parameters.

Next, we apply multi-objective optimization to characterize the Pareto frontier for combinations of three or more design parameters, by following the scalarization and bi-level decomposition approach described earlier. The Pareto frontier obtained by adding the design parameter C_p to the pair (R_f,L_f) is displayed on the bottom-left plot in Fig. 9. Also shown on this plot are the results of the Sobol’ sampling used to cross-check the Pareto points determined by the Nelder-Mead simplex algorithm.

A first observation is that the performance of the NIFTE is clearly enhanced by the addition of an extra design parameter to the optimization problem, especially in terms of the maximal power output which improves by a further 15%. Values of the triplet $\left (R_{\mathrm {f}}^{+},L_{\mathrm {f}}^{+},C_{\mathrm {p}}^{+}\right)$ at points A–E along the frontier (embedded table in bottom-left plot) suggest a radically different design strategy compared with the (R_f,L_f) case on the top-right plot. Here, reaching the Pareto frontier entails decreasing the value of $R_{\mathrm {f}}^{+}$ toward zero and moving along the frontier from A to E by reducing $L_{\mathrm {f}}^{+}$ from 0.34 to 0.28 and $C_{\mathrm {p}}^{+}$ from 4.7 to 3.9 simultaneously. Notice, in particular, the high sensitivity to changes in L_f and C_p, whereby the full frontier from A to E is described by a mere 20% variation in their values. A practical implication of this high sensitivity is that mathematical models capable of highly accurate predictions will likely be paramount to support the detailed design of NIFTE, also calling for dedicated experimental campaigns to further validate these models.

It is important to note, at this point, that the selection of alternative triplets as decision variables in the multi-objective optimization problem, including (R_f,L_f,L_p) and (R_f,L_f,C_d), produces similar improvements of the NIFTE performance. This, in turn, motivates the consideration of larger parameter subsets in the multi-objective optimization. The bottom-right plot in Fig. 9, for instance, shows the Pareto frontier corresponding to the quadruplet (R_f,L_f,C_p,C_d). Comparison with the bottom-left plot shows marginal improvement in terms of exergetic efficiency or power output after adding the parameter C_d. The consequence of C_d being redundant with the other parameters (R_f,L_f,C_p) also manifests through the lack of consistent trends in the design parameter values at points A–E along the frontier (embedded table in bottomright plot). The consideration of alternative quadruplets to optimize does not produce further improvement for this NIFTE prototype. And even the optimization of (R_f,L_f,C_p,L_p,C_d,L_d,K_th) all together leads to either a maximal efficiency close to 1.2% or a maximal power output of 0.35 W. Overall, these results establish that near-optimal performance of the NIFTE can be achieved through targeted modifications of three parameters only, such as (R_f,L_f,C_p). In other words, optimizing the length and/or diameter of the feedback and power cylinders would be sufficient in practice for NIFTE to reach its full potential.

Conclusions

The NIFTE is a simple and robust technology for low-grade heat conversion to useful power. Early-stage prototypes were found to exhibit relatively low thermal and exergetic efficiencies, but these two-phase thermofluidic engines are appealing compared to other low-grade heat recovery technologies since they can operate with temperature differences as small as (and possibly lower than) 30 ^∘C between heat source and sink, and may be technoeconomically viable in applications where conventional heat-recovery technologies are not feasible. The main objective of this paper was, therefore, to conduct a model-based assessment of the NIFTE technology under cyclic steady-state operation, with a view to better understanding its operation and quantifying its true potential.

A nonlinear dynamic model of the NIFTE describing cyclic steady states as stable limit cycles was used as the basis for our investigations, and we extended it to account for irreversible thermal losses in nominally adiabatic components of the device. For the first time, we could predict the co-existence of multiple limit cycles for certain NIFTE configurations, which may be either stable or unstable.

A parametric analysis of the most sensitive design parameters of the NIFTE device (feedback, displacer and power cylinders, and heat-exchangers) was conducted, followed by a systematic multi-objective optimization to characterize the Pareto frontier of maximal exergetic efficiency versus maximal power output. Our results establish that failure to account for irreversible thermal losses can grossly overpredict the exergetic efficiency and power output that may be delivered by a NIFTE. The assessment conducted with the extended model predicts improvement greater than 50% in either the work output or the exergetic efficiency via targeted design modifications to an early-stage NIFTE prototype. The optimization study furthermore revealed that near-optimal performance of the NIFTE can be realized by adjusting combinations of no more than three parameters, such as the dimensions of the feedback tube and power cylinder.

These improvements are particularly relevant for the development of NIFTE pumps to support small-holder irrigation in the developing world [38]. Though for the NIFTE technology to be adopted for low-grade heat recovery at larger scale, minimization of the irreversible thermal losses via further design modifications would be worth investigating, given that these losses entail a dramatic reduction in the exergetic efficiency and work output. Other recommended follow-ups to our paper include the execution of dedicated experimental campaigns to further calibrate the NIFTE models, and the implementation of optimized NIFTE designs to verify the model-based predictions. In doing so, one should keep in mind the simplicity of the electrical analogy models used relative to the complexity of actual NIFTE devices. As more experimental data of NIFTE become available, it would be interesting to attempt the development of hybrid semi-parametric models that enable more reliable predictions by combining a physics-based models with data-driven components [39].