A single-reciprocating-piston two-phase thermofluidic prime-mover

We explore theoretically a thermodynamic heat-engine concept that has the potential of attaining a high efficiency and power density relative to competing solutions, while having a simple construction with few moving parts and dynamic seals, allowing low capital and operating costs, and long lifetimes. Specifically, an unsteady heat-engine device within which a working fluid undergoes a power cycle featuring phase-change, termed the ‘Evaporative Reciprocating-Piston Engine’ (EPRE) is considered as a potential prime mover for use in combined heat and power (CHP) applications. Based on thermal/fluid-electrical analogies, a theoretical ERPE device is conceptualized initially in the electrical-analogy domain as a linearized, closed-loop active electronic circuit model. The circuit-model representation is designed to potentially exhibit high efficiencies compared to similar, existing two-phase unsteady heat engines. From the simplified circuit model in the electrical domain, and using the thermal/fluid-electrical analogies, one possible configuration of a corresponding physical ERPE device is derived, based on an early prototype of a device currently under development that exhibits some similarities with the ERPE, and used as a physical manifestation of the proposed concept. The corresponding physical ERPE device relies on the alternating phase change of a suitable working-fluid (here, water) to drive a reciprocating displacement of a single vertical piston and to produce sustained oscillations of thermodynamic properties within an enclosed space. Four performance indicators are considered: the operational frequency, the power output, the exergy efficiency, and the heat input/temperature difference imposed externally on the device's heat exchangers that is necessary to sustain oscillations. The effects of liquid inertia, viscous drag, hydrostatic pressure, vapour compressibility and two-phase heat transfer in the various engine components/compartments are examined, via changes to thermodynamic/thermophysical/transport properties and also geometrical features of the ERPE. It is found that for high efficiency and power output: (1) the vapour dead-spaces must be minimized; (2) the length of the tube that connects the displacer and working cylinders must be of significant length; and, (3) the heat-exchanger blocks must have a low thermal resistance and high heat capacity. The methodological approach implemented in this study can be used to guide the proposal, early-stage design and verification of these complex unsteady thermodynamic systems, while offering important suggestions for improved performance and system optimization.

exhibits some similarities with the ERPE, and used as a physical manifestation of the proposed concept. The corresponding physical ERPE device relies on the alternating phase change of a suitable working-fluid (here, water) to drive a reciprocating displacement of a single vertical piston and to produce sustained oscillations of thermodynamic properties within an enclosed space. Four performance indicators are considered: the operational frequency, the power output, the exergy efficiency, and the heat input/temperature difference imposed externally on the device's heat exchangers that is necessary to sustain oscillations. The effects of liquid inertia, viscous drag, hydrostatic pressure, vapour compressibility and two-phase heat transfer in the various engine components/compartments are examined, via changes to thermodynamic/thermophysical/transport properties and also geometrical features of the ERPE. It is found that for high efficiency and power output: (1) the vapour dead-spaces must be minimized; (2) the length of the tube that connects the displacer and working cylinders must be of significant length; and, (3) the heat-exchanger blocks must have a low thermal resistance and high heat capacity. The methodological approach implemented in this study can be used to guide the proposal, early-stage design and verification of these complex unsteady thermodynamic systems, while offering important suggestions for improved performance and system optimization.

Motivation
A variety of drivers have been promoting an enhanced interest in the decentralized combined production or co-generation of heat and power (CHP). These drivers include a desire for the secure supply of sustainable primary energy, health and environmental concerns arising from the emission of combustion gases to the 5 atmosphere, the increasing penetration of intermittent renewables into the electricity grid and the aspiration to provide heating and/or power to remote areas with no or intermittent energy supplies [1,2]. The EU has been addressing the issue of energy security, and one of the priorities has been to improve energy efficiency through the creation of a decentralized network of low-power cogeneration plants [3]. The decentralized framework arises from an incentive to minimize transmission and distribution losses on the one hand; and 10 on the other, it allows a greater flexibility and variation of the delivered power depending on the local demand. Conventional external-combustion prime-mover technologies, such as Rankine/steam and Stirlingcycle engines, can be relatively complex, operationally inflexible (e.g., at part loads) and expensive at small scales and low power outputs, which is why they have been largely displaced by internal combustion engines in relevant applications. Beyond issues with vibration and noise, however, conventional internal combustion 15 engines cannot directly utilize external heat sources and are designed for the direct consumption of fossil fuels.
In addition, they have relatively high maintenance requirements owing to complex mechanical arrangements and moving seals. This motivates the proposal and development of a prime mover with the following desirable characteristics: (1) efficient power generation at small scales and low power outputs between 1 and 100 kW; and (2) compatibility with a diverse variety of heat sources such as waste heat, solar, geothermal, fossil fuels, 20 etc. Relevant heat engine technologies (e.g., organic Rankine cycle, Kalina or Stirling cycles) with power outputs below 100 kW have investment costs in the order of 2−4 £/W el [1,2,[4][5][6] whereby the costs per unit power output tend to increase at very low power outputs. Therefore, at these low power ranges, the prime mover must be extremely cost-efficient for the CHP system to be economically viable.
Thermofluidic oscillator (TFO) engines (or converters) promise to be a particularly cost-efficient solution 25 to this challenge. TFOs are a class of unsteady heat-engine within which persistent periodic thermodynamic oscillations are induced from static (steady) external thermal conditions [7]. In practice, a TFO is a network of interconnected tubes and chambers that features no or few mechanical moving-parts and seals. Thanks to its simple construction and operation, a TFO typically has low capital costs and low maintenance requirements.
The working fluid is subjected to time-varying heat transfer by alternating thermal contact with a pair 30 of hot and cold heat exchangers. This alternating heat-exchange process induces thermodynamic property oscillations of pressure, volume and temperature, giving rise to an oscillatory fluid motion. Consequently, an unsteady thermodynamic cycle is sustained that is capable of converting thermal energy to hydraulic or pumping work with the use of a (liquid or solid) piston [8]. A variety of heat sources and sinks can be used to provide and extract heat to and from the heat exchangers. TFOs can take the form of single-/gas-phase 35 systems such as thermoacoustic [9][10][11][12][13][14] and Fluidyne engines [15][16][17][18], or two-/vapour-phase systems such as the 'Non-Inertive-Feedback Thermouidic Engine' (NIFTE) [7,[19][20][21][22][23][24][25].

Objectives
In this paper we conceptualize and investigate theoretically a two-phase TFO that we term the 'Evaporative Reciprocating-Piston Engine' (ERPE). In two-phase TFOs, such as the ERPE and the aforementioned 40 NIFTE, periodic heat addition to and rejection from a working fluid lead to alternating phase-change (evaporation and condensation), and consequent pressure and volumetric displacement variations within the device.
Phase change allows TFOs to operate across low temperature differences between the heat source and sink, and facilitates small heat exchanger surfaces per unit heat transferred. During the heat addition (evaporation) phase, the mass of vapour inside the device increases leading to an increase in the pressure of the 45 vapour volume. This, in turn, leads to a positive displacement of a liquid piston into a hydraulic load.
Conversely, during the heat rejection phase, condensation in the vapour space results in a pressure reduction which produces a suction stroke. These alternating pressure and volume oscillations produce hydraulic work that can be harnessed to drive a generator at the load.
The goal of this research is to develop a theoretical ERPE concept which has the potential to exhibit a 50 higher efficiency and frequency compared to competing TFOs. The approach follows the methods employed earlier by Backhaus and Swift [12,13], Ceperley [26], Smith [7,27], Markides et al. [19,23,24] and Solanki et al. [20][21][22] according to which analogies are drawn between thermofluid processes in the device and passive electrical components. This method was validated successfully in the context of two-phase TFOs against experimental data obtained from a NIFTE prototype [20][21][22], and has proven to be useful in promoting an 55 improved understanding of the operation and performance, and also the early-stage design and optimization [25], of this technology. Initially and based on the electrical analogies, the ERPE concept is defined in the analogy domain as an electronic circuit model representation of an engine of which a physical description is still undefined. Subsequently and based on an existing device that bares some degree of similarity to the ERPE in the electrical analogy domain, approximate values of important electrical parameters and a 60 possible physical configuration are taken as indications of the design space within which an ERPE may be realised. Consequently, a physical design of a possible ERPE oscillator is derived from the electronic circuit and used in the remainder of the study. This approach is a reversal of the conventional approach whereby an existing physical representation of an engine (e.g., an actual prototype or an existing design) is modelled to understand its behaviour and predict its performance. By reversing this process and starting in the electrical 65 analogy domain, we can from the very beginning design an engine which fulfils specific desired criteria-in our case, a comparably high efficiency and power density.
To define the values of the components in the ERPE circuit and thereby develop a physical ERPE device, the geometric parameters and thermodynamic/thermophysical fluid properties relating to a TFO engine prototype currently being tested at the Boreskov Institute of Catalysis (BIC) are adopted to the new concept 70 [28]. This prototype, which is an evolution of one originally developed by Encontech B.V. (www.encontech.nl) [29], has been selected as it bears several similarities to one of the possible physical representations of the ERPE, as proposed here. Additionally, we use experimental data of pressure variations provided by BIC of their prototype as indicative values to assess if the present model predictions are realistic. Two different sub-models of the heat transfer processes at the heat exchangers of the device are compared. The first model 75 assumes a linear temperature profile (LTP) along the height of the hot and cold heat exchanger surfaces.
The second model accounts for the dynamic ability of the heat exchanger walls (DHX) to store and release heat periodically by considering explicitly the unsteady energy balance that describes the heat addition (or rejection) from the external heat source (or sink) [22]. Subsequently, a parametric study is performed in which different component characteristics are varied to examine their influence on the engine's performance 80 indicators, i.e., the operational/oscillation frequency, power output, efficiency, and necessary heat input for sustained operation/oscillation. which have yet to be defined. A possible physical description of the ERPE is shown in Fig. 8 which has been derived to be as close as possible to the BIC prototype [28]. This physical representation will be used for the whole scope of the study. The geometry and thermodynamic properties of the BIC prototype will be imposed on the physical representation of the ERPE where possible to provide realistic values for the 95 electrical elements. In the interest of comprehension, the electronic circuit model is described in the following sections as being derived from the physical representation, when in actuality, the physical representation of the engine is derived from the circuit. Using the physical representation as a starting point is the more conventional approach, and it is therefore the easier approach for the reader to understand.

ERPE Configuration and Operation
The volume at the top of the displacer cylinder in Fig. 8 contains working fluid in the vapour phase, 100 which behaves as a gas spring. This section is surrounded by the hot heat exchanger (HHX) which adds heat and evaporates the working fluid-in this case water. Below the HHX is the cold heat exchanger (CHX) section, which extracts heat from the working fluid and condenses it. Connected to the displacer piston is a mechanical spring which moves freely with the piston. The working cylinder, on the left-hand side, consists of a gas spring (argon gas) entrapped by a liquid column. The two cylinders are connected by a 105 tube equipped with an adjustable valve through which the engine's power output can be assessed. When the valve is completely open, the engine is considered to run without a load.
Since the ERPE is an oscillator, understanding the processes undergone in one reciprocating cycle fully explains its operation. Assuming a cycle to start at the top dead center (TDC) of the solid piston, the liquid-vapour interface is in contact with the HHX such that heat is added to the working fluid. The 110 working fluid evaporates, generating vapour and leading to a pressure increase in the vapour gas spring.
This pressure drives a downward acceleration of the liquid level and of the displacer piston which overshoot their equilibrium positions (halfway between the two heat exchangers) due to their inertia. The vapourliquid interface and the displacer piston always move in the same direction albeit not necessarily at the same speeds. The displacer piston is assumed to be always completely immersed in water. Water flows through 115 the load in the connection tube-thereby dissipating work-into the working cylinder elevating its water level and compressing the argon gas spring. A hydrostatic pressure difference is created between the water levels of the two cylinders. Gradually, the CHX surfaces are exposed to the hot vapour, which begins to condense. The condensation process leads to a decrease in pressure in the vapour region. The increasing hydrostatic pressure difference between the displacer and working cylinders, the increasing compression of 120 the mechanical and gas springs and the decreasing pressure in the vapour region exert restoring or suction forces on the piston and the water column. However, for a limited time, the piston and water column continue their downward displacement due to their momentum, while gradually decelerating.
As the displacer piston reaches bottom dead centre (BDC), the restoring hydrostatic, spring and gas compressibility forces gradually lead to a reversal of the piston and flow direction. The water level and 125 the displacer piston accelerate upwards and overshoot their equilibrium positions due to their inertia. Once again, the hydrostatic pressure difference and the gas and mechanical springs create restoring forces; this time, in the opposite direction, thus gradually decelerating the rising piston and water level until they reach TDC where their movement is reversed. In this manner, a full cycle is completed.

130
The development of the ERPE electronic circuit model will be described in detail below. It will be presented as a derivation from a physical representation of the ERPE in order to make it easier for the reader. In this study, however, the starting point for the conceptual development of the ERPE was the electrical-analogy model that was transformed into a physical device at later stage, as mentioned in the previous section.
If a physical engine representation exists, a set of linearized, spatially lumped sub-models can be derived Specifically, thermal resistance, fluid drag and viscosity are accounted for by introducing resistors R, gravity and gas compressibility are described by capacitors C, and fluid/piston inertia by inductors L. The voltage across an element is equivalent to a pressure P at the corresponding engine location, whereas the current is equivalent to the resulting volumetric flow-rate U through it. At the heat exchangers, the voltage corresponds to a temperature and the current corresponds to an entropy flow-rate; both are converted into an equivalent pressure and volumetric flow-rate as explained in previous work [19 -24]. Applying Gauss's, Faraday's and Ohm's laws in the electrical analogy leads to the following three fundamental equations between pressure P and flow-rate U which, with the help of the principle of superposition, can account for more than one process or effect within a particular component: or in the Laplace domain: where s = σ + iω is the Laplace variable. are small fluctuations around their respective time-mean values, e.g., P P . Therefore, for simplicity, the prime symbol will be henceforth dropped, and any variables (pressure, temperature, flow-rate, etc.) mentioned will refer to the fluctuating, time-varying components of these quantities.

155
The definitions and nominal values of the electrical components are listed in Table 1. Their derivations can be found in Refs. [19][20][21][22]. To calculate the nominal values, the geometric dimensions of the BIC prototype are adopted along with relevant properties of the presently selected working fluid, water.

Phase-change/thermal domain and heat exchangers
The liquid columns in the cylinders or tubes can be described with a resistance representing viscous drag, an inductance for the column's inertia and a capacitance for the hydrostatic pressure difference (not for the connection tube as it is horizontal) which can be deduced from a force balance on the liquid column: where µ and ρ wf are the dynamic viscosity and density of the liquid working fluid, l 0 is the equilibrium length of the liquid column, d is the diameter of the cross-section, A is the cross-section area and g is the gravitational acceleration. The Reynolds and Womersley numbers of the flow are sufficiently low to assume quasi-steady, viscous laminar and fully developed flow. Flow losses at sudden contractions and expansions are neglected. Gas springs are considered adiabatic and reversible (isentropic), and are modelled by capacitances that can be derived by linearizing the isentropic perfect-gas relations. The resulting capacitance is: Here, V 0 and P 0 are the gas spring equilibrium volume and pressure, and γ is the ratio of specific heats of the gas or vapour. The heat exchanger components (see Figure 8) stem from a heat balance at the heat exchangers, the change in entropy due to phase change and the Clausius-Clapeyron relation. This leads to a resistance and capacitance (in the DHX model; see Section 3.6 and Refs. [21,22]) in series: where ρ v,0 is the density of the vapour at equilibrium state, T 0 is the equilibrium temperature, ∆s fg is the 160 entropy of vaporization, h is the heat transfer coefficient, A s is the active surface area over which phase change heat transfer occurs, (dT /dP ) sat is the change of saturation temperature with pressure in the twophase region, m hx is the mass of a fixed part of the heat exchanger wall that participates in the heat transfer process and c hx is the heat capacity of the heat exchanger walls. The definitions and nominal values of all electrical components are listed in Table 1. Their derivation can be found in Refs. [19][20][21][22]. 165

Piston and leakage
The piston can be considered as comprising two parts. In the lower part, the piston is guided by the slide  1.23 × 10 5 1st leakage flow resistance R l,1 = 128c 2 hpµ /πc 1 c 3 6.29 × 10 7 2nd leakage flow resistance R l,2 = 128c 2 hpµ /πc 1 (c 1 − 2c 2 dp 2 ) 1.37 × 10 6 Slide bearing piston resistance 7.11 × 10 7 Piston inertia Lp = 32c 2 mp /π 2 dp 2 c 1 7.44 × 10 6 Slide bearing piston inertia these impedances can be manipulated such that the circuit consisting of the components presented in Fig. 9 and Eq. 6 are derived. A physical and intuitive understanding of these components in isolation is difficult.
For example, the division of U l in Fig. 9 into two parallel branches, U l,1 and U l,2 , cannot be understood which describe the drag that results from the movement of the piston and the linear velocity profile of the lubricating film. Eventually, one arrives at the components from Fig. 9 and defined in Eq. 6.
µ is the dynamic viscosity of the working fluid, m p is the mass of the piston, k ms is the constant of the mechanical spring, is the length of the slide bearing, ρ ss is the density of the piston (stainless steel) and δ ≈ 0.1 mm is the thickness of the lubricating film in the slide bearing.

Load model 195
To assess the performance (e.g., efficiency and power output) of the ERPE, a load must be included in the circuit model. This can be done by introducing a linear resistance at the physical location where energy can be dissipated, which can be considered as representing useful work [19][20][21][22][23][24]. However, if one wants to use data from the BIC prototype tests as indicative values for comparison and validation of the ERPE model, the model must capture the behaviour of the actual load of the BIC experimental prototype (in fact, a check-valve arrangement) as closely as possible. Pressure measurements upstream and downstream of the adjustable load valve have been performed at the BIC. These indicate that this valve cannot be accurately represented by a simple resistance. Figure 10 depicts the ratio of the amplitudes of the measured pressures downstream and upstream of the valve |P down /P up | in the frequency domain acquired with an FFT algorithm.
A single resistance would return an approximately constant amplitude ratio for all frequencies. In Fig. 10, the region in the vicinity of the oscillation frequency observed by the prototype while the measurements were performed (0.2 Hz) is of particular interest. The spectral data is found to follow the trend: which is also represented in Fig. 10. In Eq. 7, s = iω is the Laplace variable which is in this case purely imaginary as the measurements are carried out while the engine is oscillating at marginal stability. In addition, ξ, λ 1 and λ 2 are constants that are fitted such that the relation follows the spectral distribution as closely as possible in the vicinity of the oscillation frequency.
Furthermore, pressure drop measurements across the load valve, when applying steady (non-oscillating) flows, can be used to calculate a steady-flow resistance R sf value with: where R sf is the experimentally measured steady-flow resistance which depends solely on the valve setting and not on the operating frequency. It is the resistance to steady flow due to viscous drag in the valve.
Equation 8 can then be rewritten in the form: In Eq. 9, R 0 is a yet to be determined constant with the same unit as R sf . Solving Eq. 7 for P down and substituting the solution into the left hand side of Eq. 9 leads to: The expression on the left-hand side in the square brackets in Eqs. 9 and 10 is dimensionless, therefore these equations can be written in the form 1 Z P = U : which can be rewritten as: where Z lo,1 is the impedance describing the load and Z lo,2 is the impedance describing all the components downstream of the load: Finally, R 0 is evaluated from Eq. 9 for s − → 0. Thus: In the above relationships, the load impedance Z lo,1 equals the steady-flow resistance R sf if there is no 200 oscillation (s = 0). Note that P up − 0 is the pressure drop across the load valve and the whole working cylinder including the gas spring. Therefore, the sum Z lo,1 + Z lo,2 accounts for the components (resistors, inductors and capacitors) of everything downstream (to the left) of the upstream pressure sensor in Fig. 8 (mainly the working cylinder). The implementation of the load model of Eq. 12 in the electrical circuit is done by placing the two impedances in series, as is highlighted in Fig. 7(b).

Liquid column below piston
The section at the bottom of the displacer cylinder is modelled by defining a control volume of fixed size moving vertically with the solid piston. This defines a control volume that contains a time-varying mass that consists of the fluid (always in the liquid phase) under the solid piston, while accounting for the different average speeds of the flow and piston. The force balance for this section takes the form: where the subscript 'd,b' denotes values for the liquid column in the displacer cylinder below the piston.
So, e.g., m d,b is the (variable) mass of the liquid column in the control volume, a dl,b is its acceleration, is the cross-section area, S d,b is the lateral area equal to the surface area of the cylinder, du/dr| r=rc is the velocity gradient at the cylinder wall (r = r c ) and u rel is the velocity of incoming/outgoing flow (averaged across the cross-section area) relative to the velocity of the liquid column. The final term of the force balance accounts for the momentum of the mass being introduced into or out of the control volume.

Liquid column above piston
Similarly to the treatment of the section below the solid piston, the liquid column above it is modelled by defining a control volume with a time-varying mass. However, in this case, the top and bottom boundaries of the liquid column are both moving independently of each other. The upper boundary is located at the liquid-vapour interface moving with velocity u and the lower boundary is located at the top of the piston moving with velocity u p . The force balance on this region is written in a similar manner to that for the liquid column below the piston (in the previous section): The subscript 'd,a' denotes here values relating to the liquid column in the displacer cylinder above the piston, and u rel is the velocity of the incoming flow relative to the velocity of the column u. Applying the same steps as in the last section leads to the expression: Equations 16 and 18 obtained for the liquid column above and below the solid piston are added together to get the total pressure drop across the entire liquid column in the displacer cylinder (excluding the piston) which can be written in terms of a set of RLC parameters as: 3.6. Solution  Table 1. Details on these models and their solution methods can be found in Refs. [21,22].  Table 2 and R th and C hx from Table 1. Finally, the exergy efficiency can be calculated with the following formula [19-23]:

Comparison of model and experiments
The BIC prototype was tested under eight different operating conditions-three without and five with a load-which are summarized in Table 3. No temperature or heat transfer measurements were carried out along the walls of the heat exchangers. Therefore, no precise experimental values exist of the temperature gradient dT hx /dy (required by the LTP model) or the heat-input rate gradient dQ hx /dy (required by the  Table 3.  no-load pressure case in Table 3, and (b) an engine configuration featuring a load with an 18.1 bar mean pressure that matches the highest (DHX) efficiency load case in Table 3. The product of exergy efficiency and frequency η ex × f is proportional to the hydraulic power output per unit heat input, and therefore high values for this product are favourable. Similarly, low heat-input rates/temperature gradients are also favourable because they are measures of the thermal input/temperature difference required to sustain the 300 operation of the device. Each RLC parameter considered here is directly linked to the geometric, thermodynamic/thermophysical or transport properties of the ERPE that result from solid materials and/or fluid substance selection decisions and the operational characteristics of the device, as defined in Table 1.
The results of this parametric exercise can be found in Figs. 1 and 2, where only influential parameters are shown and parameters with negligible influence have been omitted. In these figures, the RLC parameters 305 are given in normalized form (e.g., R * th ) meaning that the varied parameters (e.g., R th ) are divided by their respective nominal values (e.g., R th,nom ). The nominal values are based on pre-defined geometric, thermodynamic/thermophysical and transport properties taken from the BIC prototype (see Table 1) in order to give the model a physically realisable starting point for this exercise. The considered variation range for each parameter is between 10 −3 and 10 3 times the nominal starting value.

310
Beginning with the operational frequency of a load-free ERPE shown in Fig. 1(a), the three most influential parameters are the capacitance of the argon gas spring in the working cylinder C g , the capacitance of the vapour-phase gas spring at the top of the displacer cylinder C v and the inductance of the liquid in the connection tube L c (see Fig. 8). All other parameters, including C d , C hx and R th that are also shown in Fig. 1(a), have very little influence on the frequency of the ERPE.

315
The inductance L c has the largest magnitude of all the inductances in the model while the two gas-spring capacitances C v and C g have the smallest magnitudes of all the capacitances. Dominant capacitances have small values because the corresponding impedance is given by Z = 1/sC. Therefore, changes to these three parameters appear to have the greatest influence on the frequency of the ERPE as a whole. To maximize the frequency, the two gas-spring capacitances must be reduced as much as possible. Reducing the capacitances 320 of the gas springs is analogous to decreasing the time-mean dead volumes in these parts of the device.
The dependency on L c also suggests the implementation of a reduced connection-tube inductance. Table 1 suggests measures that can be taken in the desired direction; to reduce this inductance, the cross-sectional area of the connection tube must be increased and/or the length of the tube must be shortened. Figure 1(b) shows the heat input rate gradient dQ hx /dy as a function of: the hydrostatic capacitance 325 in the displacer cylinder C d , the working cylinder gas-spring C g , the heat exchanger capacitance C hx , the vapour gas-spring capacitance C v , the connection-tube inductance L c and the thermal resistance R th . With the exception of C d and L c , reducing all of these parameters below the nominal design values defined by the BIC prototype by approximately one order of magnitude (but not much further), leads to an advantageous reduction in the necessary heat-input requirement. The thermal-domain parameters, heat-exchanger capac-330 itance C hx and thermal resistance R th , attain a minimum heat-input at slightly lower values. The thermal resistance R th is inversely proportional to the heat-transfer area and the heat transfer coefficient of phase change, so increasing these (e.g., through the addition of fins, etc. [30]) reduces the thermal resistance. Also, the heat-exchanger capacitance C hx depends proportionally on the heat capacity c hx of the heat exchanger material and the mass m hx of the portion of the heat exchangers which interacts thermally with the working 335 fluid. To reduce these, a heat exchanger material can be used which has a lower specific heat capacity.
The two most dominant parameters are the hydrostatic capacitance in the displacer cylinder C d with which the necessary heat input dQ hx /dy increases linearly (in this log-log plot), and the inductance of the liquid in the connection tube L c which allows low heat inputs when increased by at least a factor of ∼100 from its nominal value. According to Eq. 3, these effects can be achieved in practice by a decrease in the 340 cross-sectional area of the displacer cylinder, and therefore a decrease in the quantity of working fluid per unit height of the heat exchanger (requiring less thermal energy for evaporation per unit height), which decreases the corresponding hydrostatic capacitance (C d ). The only significant inductance L c can be increased, for example, by lengthening the connection tube or using a smaller tube diameter. Figure 1(b) also shows sharp inflection points or discontinuities in the trends of C g , C v and L c . As 345 explained in Section 3.6, the closed-loop transfer function has seven poles: one real and three complex conjugate pairs. For marginal stability, one of these pairs is purely imaginary while the others have negative real parts. At the discontinuities, the dominant pair of poles switches to a different pair which leads to a sudden change in the oscillating frequency and gain, and hence the heat transfer rate gradient. Discontinuities have also been observed in previous linearized modelling studies of the NIFTE [22,23].

350
A similar parametric study is carried out for an ERPE engine but now with a load and a mean pressure of 18.1 bar. The influence of various parameters on the performance indicators (frequency f , exergy efficiency η ex , heat-input rate gradient dQ hx /dy and the product η ex × f ) is examined. Of the 19 variable parameters for the load case (see Fig. 7(b)), only five influence the performance indicators significantly, namely the thermal resistance R th , the heat exchanger capacitance C hx , vapour compressibility in the displacer cylinder 355 C v , the hydrostatic capacitance in the displacer cylinder C d and the inductance of the connection tube L c .
All other parameters affect the performance indicators only marginally, even if they are varied by several orders of magnitude. Figure 2 summarizes the relevant results from this parametric study.
Similarly to the load-free case, the two heat exchanger parameters R th and C hx influence the frequency in a similar manner, as shown in Fig. 2(a). When reduced below their nominal values, these two parameters 360 lead to a monotonic increase in the frequency within the entire investigated parameter ranges. Furthermore, as with the load-free case, a low vapour-spring capacitance C v results in a significantly higher frequency, as does a low connection-tube inductance L c although to a smaller extent.
The exergetic efficiency in Fig. 2(b) is mainly influenced by the vapour-spring capacitance C v and, to a lesser extent, the connection-tube inductance L c . Reducing the capacitance C v leads to a substantial increase 365 in the efficiency, which is shown to reach theoretical values up to 90%, while increasing L c to approximately 100 times its nominal value increases the efficiency from 2% to 20%. The remaining three parameters discussed here only alter the efficiency by a few percentage points. Although the ERPE has low efficiencies at the nominal values derived from the BIC prototype, Fig. 2(b) shows that it is possible with relatively simple measures to substantially increase the efficiency (this can also be said of the frequency). Increasing 370 the pressure or reducing the volume of the vapour region in combination with increasing the length of the connection tube and reducing its diameter would increase the efficiency. In Fig. 2(b), the ERPE demonstrates a potential of exhibiting high exergetic efficiencies, and it is not expectd to be advsersely affected by the same degree of parasitic condensation that has been observed in the NIFTE [21]. The remaining of the five aforementioned parameters, C d , C hx and R th , have little to no effect on the exergy efficiency.

375
Figure 2(c) shows the behaviour of the necessary heat input rate gradient dQ hx /dy. The thermal resistance R th , the heat exchanger capacitance C hx and the vapour compressibility C v all show a similar trend: the necessary heat input reaches a minimum for values slightly below the employed nominal values. Once again, the most influential parameter is the hydrostatic capacitance in the displacer cylinder C d with which the necessary heat input dQ hx /dy increases as a power law (approximately linearly in the log-log plot) for 380 the entire range of investigated values around the nominal BIC design. Furthermore, dQ hx /dy decreases to a minimum at an L c value of around 100 times its nominal value.
The final plot in Figure 2(d) shows the behaviour of the hydraulic power output per unit heat input.
Changing the inductance L c has a comparatively small effect on the normalized power output while changes in C d have almost no effect. The normalized power output is maximized for values of the thermal resistance 385 R th a little below the nominal value (0.25×R th,nom ) and decreases for values lower and higher than that. The normalized power output also increases monotonically with the capacitance C hx before levelling off around the nominal value. The most dominant parameter is the gas spring capacitance C v with which the power output decreases sharply. This parameter should therefore be kept as small as possible.

Thermodynamic cycle diagrams
Figures 3 shows pressure-volume (P lo,1 -V lo,1 ) diagrams at (across) the load, which relate to the useful power delivered to the load; and Fig. 4 shows temperature-entropy (T -S) diagrams at the heat exchangers of the engine, which relate to the exergy input to the device and to the working fluid, and help explain the trends of the exergy efficiency and normalized power output in Fig. 2(b) and (d). The areas enclosed by the P -V diagrams (Fig. 3) represent the net work output in one oscillation/cycle that can be multiplied with the 395 frequency to obtain the hydraulic power. The T -S diagrams (Fig. 4) show the temperatures of the solid heat exchanger surfaces that the fluid contacts T hx and of the working fluid T v contained within the displacer cylinder that takes part in the thermal interaction with the heat exchangers. The areas enclosed by T hx (solid or dash-dotted lines) correspond to the exergy input to the engine, and the areas enclosed by T v (dashed or dotted lines) are the net exergy gained by the working fluid. Hence, the difference between the two areas is 400 equivalent to the exergy destruction due to heat transfer to/from the working fluid during one oscillation.
For all cases, the exergy input rate to the engine was kept constant and the time-averaged pressure was 18.1 bar. The five parameters from Fig. 2 are perturbed to 1 /5 and 5 times their respective nominal values, however, only those parameters which display significant variation in their respective diagrams are shown.
The P -V diagrams are in accordance with the power output trends from Fig. 2(d). Just like the power 405 output per unit heat input, the work output per cycle increases with increasing C hx and L c , and with decreasing R th and C v . These Lissajous plots provide additional information as to what the increase in work and power output can be attributed to. An increase in area of these curves can be caused by an increase of the amplitudes of P and V , or an improvement in the phase difference between them towards an ideal of 90 • . Figure 3(b), for example, shows us that when C v is reduced, the area (and work output) 410 increases due to a slight improvement of the phase difference and, more importantly, due to a significant elevation of the pressure amplitude (and not as much the displacement amplitude). For the three other parameters, both variable amplitudes (pressure and displacement) increase similarly and contribute to the overall improvement, especially for L c . For this latter parameter, the area increases to a greater extent than the other parameters despite the phase difference also changing substantially (Fig. 3(c)); the phase 415 differences are not strongly influenced by the variations to the thermal domain parameters, C hx and R th .
For the other two parameters, the ellipse symmetry axes both steepen at lower values. Figure 4 shows how varying C v and L c affects the exergy input, and thereby, the exergy efficiency. At high C v values ( Fig. 4(b)), the area enclosed by the Lissajous T hx -S curve (which represents the external exergy input to the device) is large compared to the curve obtained with smaller values of C v , primarily due to the 420 larger amplitudes of the entropy flow-rate. On the other hand, the area enclosed by the T v curve (the net exergy transferred to the working fluid) is virtually zero due to an extremely small temperature amplitude.
When moving to smaller values of C v (Fig. 4(a)), the phase difference between T v and S increases from 16 • for 5 × C v,nom to 73 • for 1/5 × C v,nom and the amplitude of T v increases as well. Simultaneously, the area enclosed by the T hx curve drops considerably thanks to a sharp drop in the entropy amplitude (although the 425 temperature increases here also). In summary, when reducing C v the two areas exhibit opposite behaviours: the T hx area decreases and the T v increases meaning that total exergy input to the system is smaller, while at the same the net exergy transferred to the working fluid increases. This sharp reduction in exergy destruction at the heat exchanger is the reason for the sharp increase of the exergy efficiency at lower C v .
A higher L c is also associated with an improvement in efficiency although the area of the T hx curve 430 enlarges. From 1/5 × L c,nom to 5 × L c,nom , the phase difference between T v and S rises from 26 • to 40 • , while the phase difference between T hx and S decreases slightly from 90 • to 85 • . Additionally, the temperature and entropy amplitudes associated with the T v curve increase, while for the T hx curve only the entropy amplitude increases. Thus, the area enclosed by the T v curve increases to a greater extent than that enclosed by the T hx curve. Therefore, an improvement of the exergy efficiency is observed, however it is not as substantial 435 as the increase made possible through variations of C v , because the exergy input to sustain oscillations also increases. All other parameters show no great variations in the T -S diagrams.
In summary, an efficient and high-power engine requires: (a) a very low C v , which is also the most influential parameter, (b) an R th slightly below the nominal value, (c) a nominal C hx , and (d) a moderately high L c . Requirement (a) can be implemented with a small vapour space or a high mean pressure. This also 440 leads to a high frequency. Furthermore, R th and C hx can be increased/decreased by improving/reducing the heat transfer performance or the heat transfer coefficient, whereas L c is proportional to the connection tube length and inversely proportional to the tube diameter. This configuration also leads to a high frequency, however if there is a desire to increase it further, then R th should be reduced instead of C hx if both are at their nominal values. A decrease of C hx in this case, will increase frequency but also reduce power. To 445 keep the necessary heat input low, L c should be moderately high, C d should be as small as possible (small cross-section area) and all other parameters minimize the heat input slightly below their nominal values. Parametric studies have shown that for the load-free case, the argon gas spring compressibility C g , vapour 475 gas spring compressibility C v and connection tube inductance L c -if kept at a minimum-have the highest effect on increasing the frequency. Additionally, to minimize the heat input rate gradient, the thermal resistance R th , heat exchanger capacitance C hx , displacer cylinder hydrostatic capacitance C d and vapour

Conclusions
gas spring compressibility C v should be low, while the connection tube inductance L c should be high.
In an engine equipped with a load, only five parameters prove to affect the frequency, exergy efficiency, 480 heat-input rate gradient or power output significantly: the thermal resistance R th , the heat exchanger capacitance C hx , the vapour compressibility C v , the connection tube inductance L c and the hydrostatic capacitance C d of the displacer cylinder. The most influential parameter is C v . To increase the power output, exergy efficiency and frequency, the vapour compressibility C v must be as low as possible while the connection tube inductance L c should be kept moderately high. This can be achieved by decreasing the 485 time-averaged vapour space, increasing the mean pressure, reducing the diameter of the connection tube or increasing its length. Additionally, keeping the thermal resistance R th and heat exchanger capacitance C hx at intermediate (nominal) levels will further increase power output and efficiency. The resistance is inversely proportional to the heat transfer coefficient and heat transfer area, while the capacitance is proportional to the heat capacity of the walls and their mass. Finally, to keep the necessary heat input low, the parameters 490 are optimal around their nominal values except for L c , which should be moderately high (as before for the efficiency), and the hydrostatic capacitance in the displacer cylinder C d , which should be as small as possible by decreasing the cross-section area of the cylinder.
The deterioration in efficiency and power output per unit heat input depending on C v and L c can be ascribed to an increase in required exergy input to the engine which does not translate into an increase in 495 work output. The exergy destruction mainly occurs at the transfer of heat to the working fluid, often due to an unfavourable phase difference. The insight from these results offers important suggestions and guidance for improving performance and for system optimization.

Acknowledgement
This research was performed under the UNIHEAT project. The authors wish to acknowledge the Skolkovo 500 Foundation and BP for financial support. The authors would also like to thank Mr. Christoph Kirmse from the CEP Laboratory for his invaluable contribution and input.

Appendix A: Piston and Leakage Model
To derive first-order ordinary differential equations for the leakage flow and the piston motion, the Navier-Stokes equations for the flow and a force balance (Newton's second law) on the piston are applied. Figure 5 shows a schematic of the piston and leakage model. A force balance on the piston gives: where m p is the mass and A p the cross-section area of the piston, A shell is the cylindrical surface on the vertical sides the piston in contact with the liquid, and z is the vertical displacement of the piston.

505
Since only the oscillating (fluctuating) component of the motion is considered, z = z − z, m p g drops out of the analysis (see relevant discussion on this point in Section 3): which is an equation only for the oscillatory motion but with the primes dropped as before in Section 3. Figure 1: Schematic of piston and leakage model. The piston height is denoted by hp, while dp and dc are the piston and cylinder diameter respectively. The pressure at the top of the piston is P 1 and P 2 at the bottom. The leakage flow is referred to as U l while Up is the flow-rate of the piston, and k is the spring constant.
from which it is clear that the impedance Z l , in the electronic circuit representation, can be considered as consisting of an inductance L l , resistance R l,2 , capacitance C l in series and a resistance R l connected in parallel to the former three as in Fig. 6. The details of the four components are given in Eq. 6. There is no physical interpretation for the division of the U l flow into two parallel branches, or for the associated components, as there is for the rest of the engine. These interconnections and components are a result of algebraic manipulations of the coupled Navier-Stokes and force-balance equations above.
In the lower part of the piston (subscript 'b'), the leakage flow U l is diverted into two parallel and identical 515 channels. The fluid in these channels behaves like a liquid column and, therefore, can be modelled by using an inductance to account for the liquid's inertia L b,l and a resistance R b,l to represent viscous drag. These components must be placed in the U l branch of the circuit in Fig. 6, to the right of what is shown in this figure (see Fig. 9 for the complete set of connectinos), and since there are two identical channels, there are two identical parallel such branches in the circuit. Similarly, the solid-piston movement in the slide bearing 520 is modelled by using a resistance R b,p in series with an inductance L b,p to account for the friction of the piston and its inertia. These two components are placed in piston branch of the circuit U p . Finally, one arrives at the components in Eq. 6 and the circuit representation in Fig. 9.
connectionqtubeq('c') liquidq columnq('d'),q pistonq ('p')q andq leakageq ('l')q inq displacerqcylinder argonqgasqspringq('g') (noqload) workingqcylinder liquidqcolumnq('w') (noqload) heatqexchangers/ thermalqdomainq('th') vapourqgasqspringq('v')          T hx ; nom T v ; nom T hx ; 5L c,nom T v ; 5L c,nom (d) Figure 9: Linearized T -S cycle diagrams for 5 times and one fifth of nominal values of Cv ((a) and (b)); and Lc ((c) and (d)). The areas enclosed by the solid and dash-dotted lines (T hx ) represent the net exergy made available to the device in one oscillation; the dashed lines and dotted lines (Tv) represent the net exergy gained by the working fluid in one oscillation. The difference between both cycles amounts to the exergy destruction due to irreversible heat transfer between the heat exchangers and the working fluid. The temperature values on the vertical axes are not absolute temperatures but temperature differences to the respective equilibrium temperatures.