Simulation, manufacture and experimental validation of a novel single-acting free-piston Stirling engine electric generator

Applied Energy 263 (2020) 114585

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Simulation, manufacture and experimental validation of a novel singleacting free-piston Stirling engine electric generator

T

⁎

B.J.G. de la Bata, , R.T. Dobsona, T.M. Harmsa, A.J. Bellb a b

Department of Mechanical and Mechatronic Engineering, University of Stellenbosch Private Bag X1, Matieland 7602, South Africa Sasol Energy Technology, Cape Town, South Africa

H I GH L IG H T S

third-order model of a free-piston Stirling engine is derived. • AA rigorous element analysis of the linear generator is performed. • Thefinite of a novel single-acting free-piston Stirling engine is presented. • The design • validity of the proposed numerical model is shown experimentally.

A R T I C LE I N FO

A B S T R A C T

Keywords: Free-piston Stirling engine Finite element analysis Third-order modelling Validation

Owing to its high thermal efficiency, fuel flexibility, low vibration and noise, and low emissions, the free-piston Stirling engine has in recent years attracted renewed interest worldwide for uses specifically relating to microcombined heat and power generation. To aid prospective engine designers with the modelling and analysis of such engines, this paper presents the numerical simulation, manufacture and experimental validation of a freepiston Stirling engine electric generator. The paper firstly presents a 100 W electrical free-piston Stirling engine developed at Stellenbosch University. Secondly, an overview of a fully-explicit, one-dimensional numerical model is given in which both the engine thermodynamic and kinematic behaviours are solved as an initial-value problem. Thirdly, an experiential case study of passive engine operation is presented in which a peak electrical output of between 60 and 70 W was delivered. The obtained experimental data clearly validates the numerical model, with the piston stroke deviating by 2.61%; the piston-displacer phase difference by 12.42%; the workspace indicated power by 15.15%; and the average electrical output power by 23.30%. For future work it is recommended that the validated model be used to develop a more optimised and task-suited engine. The use of an active feedback controller is also recommended so that piston-casing collision can be eliminated.

1. Introduction 1.1. Background Energy is the lifeblood of our modern society; it is a vital requirement for the functioning of nearly all forms of productive means and it impacts all sectors of our economies, from best business practices through healthcare, transportation and agriculture, to education, communication and leisure, to name only a few [1]. The depletion of earth's conventional energy resources, combined with the global increase in greenhouse gas emissions, has over the last decade led to an increased urgency to develop more renewable and sustainable energy alternatives. Owing to its high thermal efficiency, fuel flexibility, low

⁎

noise and vibration, low emissions and nearly zero maintenance requirement, free-piston Stirling engines (FPSEs) have in recent years attracted renewed interest worldwide for uses specifically relating to micro-Combined Heat and Power (mCHP) generation [2,3]. First patented by the Reverend Robert Stirling in 1816 as an alternative to the steam engine, a Stirling engine is a closed-cycle external combustion heat engine that converts heat into usable mechanical (or electrical) energy via the cyclic expansion and contraction of a hermetically sealed gaseous working fluid [4]. In 1964, Professor William Beale from Ohio University brought forward a unique modification to the Stirling engine [5], by replacing the conventional crank mechanism of a Stirling engine with flexure springs that attach to the respective pistons. Since FPSEs derive heat from external sources, they have been

Corresponding author. E-mail addresses: [email protected] (B.J.G. de la Bat), [email protected] (R.T. Dobson), [email protected] (T.M. Harms).

https://doi.org/10.1016/j.apenergy.2020.114585 Received 18 October 2019; Received in revised form 26 January 2020; Accepted 27 January 2020 0306-2619/ © 2020 Elsevier Ltd. All rights reserved.

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Nomenclature

η θ λ μ ρ τ ωn

Roman symbols

A a Ax Az Cf cv cpis d dh e F f h i k KL L m N Nu P R Rg Re T t u V v x ṁ Q̇ Ẇ ẋ

area acceleration cross-sectional area surface area in z-direction skin coefficient of friction specific heat at constant volume piston viscous damping coefficient diameter hydraulic diameter electromotive force force frequency convective heat transfer coefficient electrical current thermal conductivity, spring constant minor loss coefficient length mass number of control volumes Nusselt number pressure thermal or Ohmic resistance specific gas constant Reynolds number temperature time specific internal energy volume, Voltage velocity displacement mass flow rate heat transfer rate rate of work done velocity

efficiency phase difference flux linkange dynamic viscosity density shear stress natural frequency

Superscripts

Δt t

finite time-step duration denotes current time-step

Subscripts

0 bs c d e f g h i j k m p r s w

initial backspace compression space displacer expansion space flexure spring gas heater, hot-end, hot-side i -th numbered scalar control volume j -th numbered vector control volume, interface node cooler, cold-end, cold-side mesh piston regenerator shaft (displacer) water

Acronyms and abbreviations CFD DAQ FPSE mCHP

Computational Fluid Dynamics Data Aqcuisition Computer Free piston Stirling engine micro-Combined Heat and Power

Greek symbols

ε

void fraction

[17], which envisages using an array of fission-powered Advanced Stirling Convertor (ASC) units to produce between 1 and 10 kWe of reliable, zero-maintenance power for future planetary exploration missions. In March 2018, the first round of experimental tests on the socalled Kilopower Reactor Using Stirling TechnologY (KRUSTY) was successfully completed [11]. Breakthrough technological developments of the space industry have also been applied in the renewable energy sector. However, despite the notable advantages of high thermal efficiency, quiet operation, zero maintenance requirements and fuel flexibility, widespread commercial success of such systems have been limited [18]. In 2012, for example, Infinia Corporation, a fore-runner in the development of FPSEs who were marketing their FPSE technology for remote solar power electricity generation, filed for bankruptcy [19]. At the time of writing, the only commercially available FPSE units were that of Microgen Engine Corporation [20] and Qnergy [21], despite ongoing government-funded research in the United States [22] towards developing a natural-gas-fuelled FPSE for residential micro heat and power application.

applied to a range of applications including waste heat recovery and utilisation [6], bio-mass or methane combustion [7], solar thermal application [8,9] and even nuclear thermal heating [10,11]. Under contract from NASA and the U.S. Department of Energy, Mechanical Technology Incorporated (MTI) invested considerable effort in the development of FPSEs suited for both concentrated solar power application [8] and for use in deep space missions as the primary electrical generator [12]. As early as 1977, MTI had demonstrated the operation of a 2 kW FPSE with linear alternator [13] and by September 1993, had succeeded in developing a 25 kWe prototype that operated under a mean charge pressure of 15 MPa using helium gas as the working fluid. Also in pursuit of this undertaking, Sunpower Incorporated developed the RE-1000 Stirling engine, documented by Schreiber and Geng [14], which to this day remains one of the most extensively studied experimental free-piston Stirling engines. From 1999 to 2006, NASA's Glenn Research Center actively supported the development of a 110 We high-efficiency Advanced Stirling Radioisotope Generator (ASRG), which was envisaged as a more efficient replacement to the conventional Radioisotope Thermoelectric Generators (RTGs) that are currently being used for deep space application [15,16]. In 2015, work on the ASRG was redirected to NASA's Kilopower Development Program 2

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

his dissertation, Goldberg derives first a linearised dynamic analysis model of a back-to-back FPSE based on state-space theory, and secondly, a fully-implicit third-order model that includes a turbulence model for the engine heat exchanger sections. Under contract for NASA’s Lewis Research Centre, Goldberg [41] later validated his implicit third-order model against experimental results of NASA’s Space Demonstrator Engine. In 1992, Yuan and Spradley [42] from Lockheed Missiles and Space Company published a study on the development of a third-order model for split-cycle Stirling refrigerators – with validation having been performed on the Lucas Stirling Refrigerator. Later in 2006, Andersen [43] published a dissertation on simulating cyclic thermodynamic processes using an implicit one-dimensional model. Relating to Stirling engines, Andersen's model has shown good correlation to experimental results from both the SM5 Stirling engine and the Twinbird free-piston cryocooler. Deetlefs [44] derived a quasi-static third-order model and verified results against that of a modified Beale B-10 engine from Sunpower Inc. Wang, Dubey, Choo and Duan [45] derived an implicit third-order model for kinematic Stirling engines and has shown good correlation to experimental results from the GPU-3 kinematic Stirling engine. Thermoacoustics modelling of free-piston Stirling engines have also been employed [46] to aid in the understanding of FPSE operational mechanics. In this approach, the authors show that by matching the acoustic impedance of the thermodynamic cycle to that of the linear alternator, the piston stroke and frequency may be determined implicitly. From the published literature it is evident that most of the models developed to date consider the FPSE as a special kind of kinematic Stirling engine. Typically, a modelling methodology is followed to first estimate the displacer and piston kinematics from linearised dynamic analysis. Secondly, either an ideal-isothermal or ideal-adiabatic cycle analysis is performed to estimate the engine thermal performance, including power output and net efficiency. Third-order models may offer a considerable improvement over second-order methods as the working fluid transport equations are considered and the intrinsic loss mechanisms are retained. However, as with second-order models, there has been a preference to studying the engine thermodynamics separately from the piston kinematics and generator electrodynamics.

1.2. Overview of modelling strategies As a consequence of this poor market penetration, the published literature on experimental test data and validated numerical models remain limited to this day [23]. In fact, most of the studies found in published literature simulate the operation of Stirling engines by adopting so-called decoupled ‘first-order' or ‘second-order' modelling approaches [24,25]. In a first-order model, an ideal-isothermal loss-free cycle analysis is performed based on the Schmidt assumption to calculate the engine output power and efficiency [26]. Second-order models start with a simplified cycle analysis to determine the basic heat input and boundary work output of the engine. Either an ideal-isothermal (Schmidt) or ideal-adiabatic (Finkelstein) cycle analysis is performed, distinguishable in the way in which heat transfer is modelled in the expansion and compression space volumes [27,28]. Once the idealised engine performance is calculated, individually identifiable power loss terms are subtracted from the estimated power output and non-idealised heat losses are added to the required heat input, to more realistically estimate the engine’s net power output and efficiency [29,30]. Relating to FPSEs in particular, the displacer and piston kinematics are traditionally simplified to that of a linearised dynamic system, requiring that the assumptions of the Schmidt ideal-isothermal analysis be inferred. The advantage of employing a linearised dynamic analysis is that the system equations may be solved analytically by means of a state space analysis or control systems theory; the drawback being that the system equations do not yield a solution that is unique for both displacer and piston stroke [26]. Often credited for first proposing the use of a linearised harmonic model, Redlich and Berchowitz [31] investigated the requirements for steady operation of FPSEs. By using the ideal-isothermal modelling approach to present the working fluid thermodynamics, and transforming the linearised equations to complex space, the authors investigated regions of stable and sustainable operation for different FPSE configurations. As an improvement over the linearised dynamic approach, non-linearised dynamic models of FPSEs have also been studied. Typically these models are based on the assumptions of an ideal-adiabatic cycle analysis, requiring that the displacer and piston equations be computed numerically [23]. Using a second-order model, Karabulut [32] showed how a non-linear dynamic analysis of the displacer and piston may be solved as an initial value problem to a hypothesised free-piston Stirling engine. Following a similar approach, Formosa [33] derived a semi-analytical model of a freepiston engine in which the thermodynamic engine cycle was linked to the displacer and piston kinematics through an iterative solver. Numerical results were verified against experimental results of NASA’s RE1000 FPSE. Strauss [34] derived a second-order model based on an ideal-adiabatic cycle and investigated direct displacer and piston control strategies. His second-order model was validated against experimental results of the GPU-3 kinematic Stirling engine. TavakolpourSaleh, Zare and Bahreman [35] derived a dynamic model for a FPSE with non-linear springs based on a perturbation method, in which the thermodynamic cycle was based on the assumptions of Schmidt. In their model, the engine dynamic modelling was coupled to empirical finitetime heat transfer relations of the engine heat exchangers, allowing for more accurate estimates the steady-state internal gas temperatures to be used. Numerical results were validated against experimental results of a modified Beale-B10 FPSE. As a significant improvement over second-order models, so-called third-order theoretical models solve the conservation equations of mass, momentum and energy over a network of finite gas control volumes numerically [36]. However, few studies have been published to date, with the work done by Israel Urieli [37,38] still remaining the most thoroughly documented third-order model available in the literature. Independent from Urieli, Schock [39] from Fairchild Corporation developed a rigorous third-order computer model. In 1987, Goldberg [40] from the University of the Witwatersrand extended on Urieli's works. In

1.3. Research aim and paper outline In spite of the vast increase in computational speed that has made solving complicated third-order models quicker, no validated numerical model exists in published literature that captures the transient, nonlinearised free-running motions of the displacer and power piston in a coupled manner to solving the working fluid thermodynamics and generator electrodynamics. At the time of writing, the computer program Sage by Gedeon Associates [47] remains the only commercially available software package for modelling Stirling cycle machines. Although Sage computes the periodic steady-steady operation, by solving conservation of mass, momentum and energy equations implicitly over a discretised spatial domain, the power piston stroke is still required as a model input parameter. To address the knowledge shortfall of a validated, comprehensive numerical model that is capable of computing the time-dependent thermodynamic, kinematic and electrodynamic operating behaviour of a FPSE in a coupled manner, this paper outlines the derivation, theoretical simulation and experimental evaluation of a free-piston Stirling engine electric generator. The computer simulation program developed throughout this study requires only that the engine heater and cooling water inlet temperature be specified. Accordingly, the time-dependent engine motion and generator load power evolve naturally as simulation outputs, therefore the model is capable of emulating the real dynamic behaviour of such an engine. As no such comprehensive and validated model is currently available in the published literature, this study presents a significant contribution to the still-limited body of literature surrounding the modelling of free-piston Stirling engine operational 3

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

operating frequency of 30 Hz, a peak electrical power output of 100 W is produced. The design of the displacer motor is similar to that of the generator, but with a key distinction being that the motor magnet train consist of a single row of north polarised magnets, thereby making for a more compact but weaker design. The motor plunger, machined from PTFE plastic, attaches to a phosphorous bronze guide piston that centres the magnet train with respect to the motor core. This guide piston attaches to the displacer flexure via an adjustable connecting collar. During final assembly, the radial angle of the displacer shaft (and displacer) relative to the power piston could be adjusted so that the displacer assembly could be set up with the least amount of mechanical friction.

mechanics. It is envisaged that the use of this validated model will aid as a powerful design tool to prospective FPSE researchers, from which ultimately, an optimised and task-suited FPSE may be developed. This paper firstly gives a descriptive overview of a 100 We (peak) engine prototype that serves as a technology demonstrator. Secondly, a discretised mathematical model of the engine is presented and the main transport equations are given. Thirdly, by inferring generator characteristic equations from a transient finite element analysis of the linear generator, the generator performance was accurately emulated in the numerical simulation program. This allows the linked thermodynamic and electrodynamic behaviour of the engine to be studied in a relatively computationally inexpensive way, whilst retaining the real generator induced currents and electromagnetic forces acting on the piston. Fourthly, a descriptive overview of the experimental setup is given before lastly, a comparison between simulated and experimental data is made for a demonstrative case study of self-sustained engine operation.

3. Theoretical modelling In this paper, only an overview of the governing system equations of the numerical model SENSE (Stirling Engine Numerical Simulation Environment) is given. Refer to the dissertation by De la Bat [50] for a detailed derivation of the system equations. The modelling methodology may briefly be outlined as follows:

2. Engine description To validate the numerical model derived throughout Section 3, a 100 We (peak) FPSE prototype was designed and manufactured. This first-generation prototype has remained completely unoptimised for the scope of this study; the engine serves primarily to evaluate the numerical model's capability of emulating its underlying operating mechanism. The geometric layout of the prototype is similar to Sunpower Inc.'s P2A FPSE (formally known as the EG-100 [48]) and Microgen’s [20] commercially available 1.0 kW FPSE, but with the addition of a linear motor being attached directly to the displacer shaft. In Table 1, the operating and main design parameters of this prototype are listed. A trial-and-error design-methodology was ensued using the numerical model to identify design-specific considerations relating to the engine geometry. This includes the determination of the displacer and piston masses, spring stiffnesses, heat exchanger sizing, etc., while maintaining the objective of manufacturing an inexpensive, laboratorybased and testable prototype using readily available engineering materials and machining processes. Fig. 1 illustrates a quarter-sectioned view of the assembled prototype. The engine casing comprises of three steel compartments that seal at the flanges using silicone-rubber gaskets. The design of safety-critical components such as the engine casing, the displacer, power piston and flexure springs were complemented by finite element stress and thermal analyses performed using Siemens NX 11. The engine pressure casing has a peak allowable pressure rating of 2.0 MPa gauge when operated with helium as the working fluid. Attached to the so-called backspace pressure hub, a safety relief valve (not visible here) ensures the engine pressure remains within this limit. Heating is supplied by a 1.0 kW ceramic heating jacket mounted around the outside of the steel heater head, whereas at the cold end, heat is rejected to a radially-finned water cooling jacket. The internal heater and cooler fins were manufactured from copper by means of electro-discharge machining (wire cutting) 135 axially-aligned channels into annular cylinders to achieve increased heat exchanger surfaces. The regenerator canister, not shown in this rendering, locate-fits between the heater and cooling fins and consists of wound-wire copper-60 mesh strips, forming an macroscopic void volume of approximately 80%. The displacer and power piston are each supported axially by respective spiral flexure springs that exhibit a linear force-displacement relationship. The design of the flexure springs was based on a twolegged variant of the Oxford Aerospace Cryocooler spring layout [49], although the final design geometry was determined by aid of a finite element analysis. The power piston is bolted directly to the generator magnet plunger that contains 32 radially magnetised arc magnets (NeFeB grade 42H), making up the power piston assembly. As seen from Fig. 1 (right), the arc magnets were arranged axially in a north-southsouth-north polarisation to make up four magnet ‘rings'. When the plunger reciprocates over the design-specified 20 mm stroke at an

• Dividing the engine workspace and engine metal walls into arrays of one-dimensional finite control volumes. • Deriving time-dependent finite difference equations by applying the

•

principles of conservation of mass, momentum and energy to the gas control volumes, and by applying the conservation of energy principle to the metal control volumes. Deriving the non-linearised equations of motion for the displacer and piston by applying Newton's second law of motion to the respective pistons. Deriving finite difference equations for the generator voltage and current. Solving the resulting finite difference equations for each gas and metal control volume as an initial value problem using a fully-explicit numerical method with linear upwind differencing. Solving the displacer and piston velocities and displacements in a fully-

Table 1 Nominal operating and design parameters of the 100 We free-piston Stirling engine prototype. Value

Unit

He 600 <50 1.80 30±0.5 18–24 18–24 100

°C °C MPa Hz mm mm W

412.08 77.8 1251.34 326.66 13.63 17.57 10.52 11.81 7.61 1970 48 10 30 30 0.526 0.952

cm2 % cm2 cm2 cm3 cm3 cm3 cm3 cm3 cm3 mm mm mm μm kg kg

Displacer spring stiffness, kd Piston spring stiffness, kp

20 019 24 610

N/m N/m

Displacer sliding friction, Fd Piston viscous friction, cp

1.97 40

N Ns/m

Operating parameters Working fluid Hot-end temperature Cold-end temperature Charge pressure (absolute) Operating frequency Nominal displacer stroke Nominal piston stroke Generator output power (peak) Main design parameters Heater-fin heat transfer area (135 fins) Regenerator void fraction Regenerator mesh heat transfer area Cooler-fin heat transfer area (135 fins) Heater channel void volume Regenerator void volume Cooler channel void volume Expansion manifold dead volume Compression manifold dead volume Backspace volume Displacer bore Displacer shaft diameter Piston bore Piston leakage gap (wall-to-wall) Displacer mass, md Piston mass, mp

4

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Fig. 1. Rendered quarter-section illustration of the assembled engine prototype (left) and isolated view of the internal engine assembly (rights) with the main components as labelled.

control volume in an explicit fashion with forward marching in time to compute each cell's following time-step mass and specific internal energy. At the interface between two consecutive control volumes, where advective transport occurs, a boundary node was defined. Boundary nodes are denoted by the subscript ‘ j ’ and depicted by ‘○’ and serve to define the location, magnitude and direction of advective fluid flow that results due to an imbalance in pressure between neighbouring cells. As illustrated by the shaded control volume in Fig. 2, vector control volumes were staggered-stacked about the boundary nodes j – over half the upstream and downstream scalar control volumes. The principle of conservation of momentum was applied to each vector control volume in order to compute the fluid momentum and velocity at the boundary node, as outlined by Patankar and Spalding [51].

explicit manner using the internal engine pressures and load forces calculated from the previous time-step. The novelty of this modelling approach is that only the heater head and coolant inlet temperatures are required as model boundary conditions. Accordingly, the time-dependent engine motion (including the displacer-piston phase difference) and the generator voltage and current evolve naturally as simulation outputs (i.e. an initial value problem). Consequently, as will be shown in Section 6, the numerical model is capable of emulating realistic engine dynamic behaviour to a remarkable degree of accuracy.

3.1. Discretisation of the working space

3.2. Governing transport equations

Fig. 2 illustrates how the engine working fluid and metal walls were discretised into a number of control volumes, having cell centres denoted by i and numbered 1 to N . The heater, regenerator and cooler sections were each subdivided into a predefined number of control volumes respectively, whereas the expansion, compression and backspace chambers were modelled as lumped Eulerian control volumes with moving boundaries [38]. As a boundary condition, the electrical heating jacket may either be defined as a constant temperature heat ̇ source Th , or the heater power Qheat may be specified and the heater wall temperature calculated as a model output. At the engine cold end the water cooling jacket was subdivided into a corresponding number of control volumes so as to model the gradual increase in cooling water temperature, from inlet temperature Tw, i to outlet temperature Tw, o . The working fluid control volumes, also referred to as scalar elements or cells, were defined as containing a finite amount of mass with uniform thermodynamic properties such as density, pressure and temperature. Referring to Fig. 2, these scalar control volumes are denoted by the subscript ‘i ’ with cell centres being depicted by ‘●’. The principles of conservation of mass and energy were applied to each scalar

The governing transport of the working fluid were derived by integrating the general differential transport equations for mass, momentum and energy described by Bird, Stewart and Lightfoot [52] over the one-dimensional array of working fluid control volumes depicted in Fig. 2. 3.2.1. Conservation of mass For the one-dimensional case in x , the continuity equation is given:

∂ρ ∂ = ρv ∂t ∂x

(1)

and in discretised form, for the general scalar control volume i , is given in integral form by:

Δmi = [ρj vj Aj − ρj + 1 vj + 1 Aj + 1 ] Δt

(2)

where the density ρj is obtained from first-order upwind differencing, the velocity vj is determined from the momentum-balance equation, and 5

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

the conservation of momentum equation is given in integral (conservative) form by:

∂ ∂ ∂ ∂τ → (ρv ) = − (ρvv ) − (P ) − s + ρg ∂t ∂x ∂x ∂x

(4)

If the force of gravity is assumed negligible, and if minor loss effects are considered for sudden expansions and contractions, the discretised form of the momentum balance equation, as applied to the generalised staggered control volume j , is given by: Δ Δt

(mj vj ) = [ρi − 1 vi − 1 |vi − 1 | A xi − 1 − ρi vi |vi | A xi ] + [Pi − 1 − Pi ] Aj − ⎡τsi − 1 ⎣

A zi − 1 2

+ τsi

A zi 2

⎤ ⎦

1

− ⎡ 2 KL ρj vj 2A xj ⎤ ⎣ ⎦

(5)

where the absolute velocity term indicates that the direction of flow is retained. In this model the vector cell mass mj is taken as the sum of half the adjacent cell masses over which the vector volumes were staggered, Ax is the cross-sectional flow area perpendicular to velocity, Az is the wetted flow area parallel to velocity and Aj is the fluid-to-fluid interface area at node j . The cell centre velocity vi is inferred from a volumetric flow interpolation between neighbouring interface (node) velocities vj . The fluid shear stress experienced by a finite cell is simplified to an empirically-related shear stresses τs , defined in terms of the cell centre velocity vi by means of an empirical steady-flow skin coefficient of friction Cf :

τsi =

1 Cf ρ vi |vi | 2 i i

(6)

Fluid friction relations for the heater and cooler channels were based on smooth annular duct flow, using the appropriate hydraulic diameters and friction factors for either turbulent or laminar flow. For the regenerator channel, Gedeon and Wood [54] relate the skin friction coefficient of a generic woven-screen regenerator mesh to a modified three-parameter Ergun equation, by:

Cf = (

3.2.3. Conservation of energy For one-dimensional flow in x , the thermal energy equation is given in conservative differential form, by:

the fluid-to-fluid interface area perpendicular to velocity is used for Aj . Gas leakage between the compression and backspace volumes occur through a small annular cavity with height h that exists between the piston wall and its bore. According to Grinnell [53], if the assumption of isothermal, fully-developed laminar flow is adopted, the gas leakage mass flow rate may be estimated analytically by:

ṁ l =

24μRg TL

(Pc 2

− Pbs

2)

(7)

For gas flowing into or out of the expansion and compression space chambers, an additional minor-loss term is included in the momentum equation (i.e. the last term of Eq. (5)). This minor loss force is given in terms of a lumped minor loss coefficient KL , assumed as KL = 1.55 for sudden expansion effects (flow into either the expansion or compression manifold), and 2.55 for sudden contraction effects [55]. No additional minor loss effects have been considered at the interfaces of the heater and cooler with the regenerator, as the empirically-related regenerator friction coefficient is already inclusive of these minor losses.

Fig. 2. Cross-sectional illustration of the engine working space, showing how the working fluid channel and engine metal sections were discretised into onedimensional cells.

πDh3

129 + 2.91Re−0.103 )/4 Re

∂q ∂ ∂v ∂ ∂ (ρu) = − (ρuv ) − −P − τxx ∂x ∂x ∂x ∂t ∂x

(8)

and in discretised form for the general scalar control volume i : Δ Δt

(mi ui ) = [ρj vj Aj uj − ρj + 1 vj + 1 Aj + 1 uj + 1] ̇ + Pi [vj Aj − vj + 1 Aj + 1 ] + Qnet + Ẇb + τ f A zi |vi|

(3)

i

in which flow is defined as being positive for leakage occurring from the compression space towards the backspace.

(9)

where the specific internal energy uj is obtained from first-order upwind ̇ represents the net energy transferred to the gas condifferencing, Qnet trol volume from the engine metal and Ẇb is the additional net rate of boundary work done by the moving piston and displacer on the gas control volumes in the expansion and compression spaces. The

3.2.2. Conservation of momentum With the assumption of one-dimensional plug flow in the case of x , 6

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

∑ Fp= (Pc − Pbs )(Ap − As ) − kp (x p − x p0)

irreversible increase in thermal heat by viscous dissipation is encapsulated by the surface friction term, where the absolute velocity |v¯i | is indicative of the increase in thermal energy regardless of flow direction. The fluid-to-surface convective heat transfer rates were modelled by relating the appropriate hydraulic diameters, surface friction factors and local Reynolds numbers to empirically-related Nusselt numbers. For the gas control volumes in the smooth-duct heater and cooler channels, the Nusselt numbers were computed from the Gnielinski equation [55], and a fin-effectiveness analogy was used to simplify numerical modelling to a uniform fin-base temperature [56], computed for each control volume. For convective heat transfer between the gas and the regenerator mesh metal, Gedeon and Wood's [54] relation for Nusselt number was used:

Nu = (1.0 + 0.99Pe 0.66) εr 1.79

− Fmp − Ff , p + mp g

Applying Newton's second law and integrating the acceleration term temporally (assuming constant acceleration during a finite time-step Δt ) yields the displacer and piston displacements as: t

∑ Fd ⎞ vdt + Δt = vdt + Δt ⎛ ⎝ md ⎠ ⎜

(16)

Further integrating the velocity profile with respect to time yields the displacer and piston displacements:

(10)

t

x dt + Δt = x dt + (Δt ) vdt +

(Δt )2 ⎛ ∑ Fd ⎞ 2 ⎝ md ⎠

x pt + Δt = x pt + (Δt ) vpt +

(Δt )2 ⎛ ∑ Fp ⎞ ⎜ ⎟ 2 ⎝ mp ⎠

⎜

⎟

(17) t

(18)

As shown by the program flow diagram in Fig. 4, if the newly determined velocities lead either to the displacer or piston colliding with the engine casing, the new time-step velocity is set equal to zero before integrating to displacement. This assumption implies a completely nonelastic collision occurring. On the other hand, in the event of the displacer colliding with the power piston, the displacer is given a velocity profile identical to that of the heavier power piston before integrating to the new time-step displacements.

(11)

3.4. Generator circuit modelling As was shown in Fig. 1, the linear generator consists of a power piston assembly that reciprocates through a soft-iron stator. The reciprocation of the magnet train induces one of three flux paths through the stator. According to Faraday's law of induction, if a change in flux linkage is presented, an opposing electromotive force (emf) will be induced in the coil that is directly proportional to the rate of change of the encapsulating flux [57]. In the case of having a coil made up of N identical windings, this induced emf e is given by:

(12)

where Rg is the specific gas constant. 3.3. Dynamic equations of motion The displacer and power piston were each treated as moving rigid bodies that are ‘sprung' to the engine casing body via the mechanical flexure springs, as seen in Fig. 2. Their governing equations of motion were derived by applying Newton's second law of motion and subsequently using Eulerian time- integration to derive relations for velocity and displacement. The forces acting on the displacer are the expansion and compression space pressure forces Pe Ad and Pc (Ad − As ) respectively, the backspace (bounce space) pressure force acting on the displacer shaft Pbs As and the flexure spring force kfd (x d − x d 0) . Although an electromagnetic force Fmd is exerted onto the displacer by the linear electric motor, in this study, the linear motor has been used only to kickstart the engine into operation. The power piston is subject to the working fluid pressure force Pc (Ap − As ) and backspace pressure force Pbs (Ap − As ) , the flexure spring force kfp (x p − x p 0) and the generator load force Fmp , having a terminal voltage of VT , g . The power piston magnetic force is comprised of both open-circuit magnetic and closed-circuit electromagnetic force constituents, as discussed in Section 4. From an experimental characterisation, the displacer and piston friction forces Ff , d and Ff , p were determined. The net body forces acting on the displacer and piston in the positive x -direction (from left to right) are given by:

e = −N

dϕ dλ =− dt dt

(19)

where λ = Nϕ is the flux linkage. When a load is attached to the coil terminals, the electrical circuit is complete and a current flows. This current interacts with the magnetic field and induces the Lorentz electromagnetic force, given by the cross product:

F = il × B

(20)

where i represents the magnitude and direction of current flowing, l the total length of the conducting coil and B the vector magnetic field (also called flux density). It is this magnetically induced force that couples the generator electrodynamic behaviour via the power piston kinematics to the engines thermodynamic behaviour. Fig. 3 depicts an electrical circuit diagram of the linear generator with its terminal lead wires VT attached to an adjustable load circuit with resistive value Rload . The internal circuit represents the generator stator, showing the induced emf as a voltage source e with a coil inductance Lcoil and internal resistance R coil . The coil inductance not only induces a drop in voltage due to its reactance, but shifts the phase of induced current and Lorentz load to lag behind that of the induced emf. This phase shift may either be reduced passively by including a tuning capacitor in series with the inductance and load, or actively, by using power electronic modulation [58]. The circuit diagram shows how, for example, two electrolytic tuning capacitors C1 and C2 may be connected in parallel with two diodes D1 and D2 in a circuit arrangement that

∑ Fd= Pe Ad − Pc (Ad − As ) − Pbs As − kd (x d − x d 0) − Fmd − Ff , d + md g

(15) t

in which c v is the specific heat at constant volume, which has been assumed to remain constant over the range of temperatures, pressures and densities considered. By employing the ideal gas law as an equation of state, the new scalar control volume pressure is calculated according to:

Pit + Δt = ρit + Δt Rg Ti t + Δt

⎟

∑ Fp ⎞ vpt + Δt = vpt + Δt ⎛⎜ ⎟ ⎝ mp ⎠

in which Pe is the local cell Peclet number and ε is the mesh porosity. Using Euler's method, Eqs. (2), (5) and (9) were integrated temporally to fully-explicit forms to yield solvable equations for the control volume mass, momentum and thermal energy at a following time-step t + Δt . By dividing the new time-step vector control volume momentum by mass, the new time-step velocity v tj + Δt at the scalar cell interface may be computed at each time step. Similarly, by dividing the new time-step scalar control volume thermal energy by its mass, the cell specific internal energy uit + Δt may be determined. Accordingly, the new time-step working fluid temperature may be computed from the change in internal energy in a fully-explicit form, by:

Ti t + Δt = Ti t + (uit + Δt − uit )/ c v

(14)

(13) 7

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

dimensional control elements. From a grid independence study it was found that using 73 control volumes in each of the heat exchanger channels yielded adequately converged results, with the net thermal-toelectrical efficiency converging to within 95% of the grid-independent solution [50]. Accordingly, lumped control volumes were used for the expansion, compression and backspace volumes so that 222 finite control volumes were simulated in total. As depicted in Fig. 4, within the program main loop, the program first considers the ‘electrodynamic part’ of the model. This relates to employing the generator characteristic functions derived throughout Section 4 to compute the induced flux linkage λ , flux density B and magnetic forces Fm acting on the power piston as functions of piston location and velocity. After computing the rate of change in flux linkage, the induced voltage, generator current, Lorentz magnetic force and output power are calculated. Subsequently, after having determined the new time-step electromagnetic forces, the ‘displacer and piston kinematic part' of the model is computed – as represented by Eqs. (15)–(18). The thermodynamic part of the model is computed in four subroutines that are representative of applying and solving the heat transfer and fluid friction relations, and the conservation equations for mass, momentum and energy. Firstly, in the heat transfer and fluid friction subroutine, the local Reynolds numbers for control volumes numbered i to N are computed and related to their locally defined skin coefficients of friction, Nusselt numbers and convective heat transfer

Fig. 3. Electrical circuit diagram of the generator and load circuit.

allows the circuit to be switched on in reverse-bias as well. In theory, when the circuit's resonant frequency matches the operating frequency of the power piston, the current phase lag is eliminated and the maximum electrical power can be extracted from the resistive load. To derive a differential equation for the load current i , the diodes were neglected for steady-state operation and the capacitors modelled as a single equivalent capacitor in series, Ct . Applying Kirchoff's voltage law yields:

L

1 di + Rt i + Q = e dt Ct

(21)

where e is the induced emf, Rt is the combined coil and load resistance, i = dQ/ dt and Q is the electrical charge stored in the capacitor. Differentiating with respect to time and substituting for the term dQ/ dt yields:

L

d 2i di 1 de + Rt + i= dt 2 dt Ct dt

(22)

which is a second-order ordinary differential equation in time. When the inductance and resistance values are known, the capacitance may be chosen to match the circuit's resonant frequency. From Eq. (22), the current may be expressed in a fully-explicit form by use of forward and central differencing schemes, as: ⎛ Δt 2 ⎡ ⎢ ⎣⎝ ⎜

it + Δt =

et + Δt − et ⎞ ⎟ − Δt ⎠

it R it + t ⎤ Δt ⎥ Ct ⎦

L + Δt / Rt

+

2Lit − Lit − Δt L + Δt / Rt

(23)

To solve Eqs. (19), (20) and (23) so-called generator characteristic functions for flux linkage, magnetic force and magnetic flux density were inferred from the results of an electromagnetic finite element analysis, discussed in Section 4, in which the required variables were characterised (mapped) against magnet train location and velocity. By employing these simplifying generator characteristic functions, the linear generator may be represented and coupled to the power piston dynamical modelling in the overarching numerical model. The coupling of the thermodynamic, electrodynamic and kinematic model constituents are further illustrated in by the program flow diagram, discussed in the following sub-section. 3.5. Numerical procedure The program flow diagram in Fig. 4 gives an overview of the principal sequence of events encompassed by the numerical model. The computer program, written in Fortran, starts by importing the model operating parameters and discretising the working fluid, metal wall, regenerator mesh and water cooling channel into arrays of one-

Fig. 4. Top-level flow diagram of the computer simulation program. 8

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

λ (x )= − 3.24 × 10−8x p5 + 1.21 × 10−8x p 4

rates using appropriate empirical correlations. Secondly, the conservation of mass is applied to each scalar control volume numbered i through N + 1 to compute the ‘new' time-step cell mass and densities. This step includes adjusting the workspace and backspace gas mass and internal energy for gas leakage. Thirdly, the gas momentum-balance equation is applied to each vector control volume j to compute the following time-step nodal velocity vj . Fourthly and lastly, the conservation of energy is applied to the scaler control volumes denoted i and the ‘new' time-step gas, metal and cooling water temperatures are computed accordingly. The new time-step temperatures and pressures are then solved from the appropriate equations of state. After computing the ‘gas thermodynamic part' of the model, a conditional function evaluates whether an engine cycle has been completed before the per-cycle performance indicators are determined. The per-cycle change in internal energy (ΔEi ) of each gas and metal cell is then investigated as a measure of dynamic steady-state convergence. To reduce simulation time, a metal-temperature acceleration scheme has been employed to accelerate the evolution of the numerical solution to reaching thermodynamic equilibrium. The accelerated scheme divides the engine metal specific heat and internal energy by an ‘acceleration factor' kaccel to increase the metal's sensitivity to changes in temperature. Then, during the first Naccel engine cycles simulated, these reduced specific heats are linearly increased until reaching a factor of unity. Additionally, at the end of each engine cycle, the backspace gas and regenerator mesh metal temperatures are adjusted accordingly as step changes with the objective being to minimise their per-cycle rates of change in temperature. With the completion of a time-step Δt , the program stores the ‘new' time-step computed variables in place of the old. The program then evaluates whether it is required to write the output variables to storage before subsequently progressing to the following time-step. When the final simulation time is reached the program is terminated.

+ 2.69 × 10−5x p3 − 3.93 × 10−6x p2 − 8.94 × 10−3x p + 2.27 × 10−4

(24)

where the flux linkage λ is given in Weber and x p the magnet displacement is given in millimetres. Fig. 6 depicts the open-circuit magnetic force exerted by magnet stators onto the moving magnet train as a function of magnet train displacement. This force significantly influences the resonating frequency of the power piston assembly, as it acts on the magnet train as a non-linear magnetic spring. Between the displacement amplitudes of −4 and +4 mm, the magnetic force exhibits a nearly linear behaviour with increasing displacement amplitude, inducing an effective spring stiffness of approximately 10.35 N/mm. Between the piston displacements of 5 mm and 10 mm, the magnetic force decreases, crossing the zero-axis magnetic force around 10 mm. It is furthermore seen that as the piston stroke is increased, the magnetic force path widens. The enclosed area within a given force-displacement path represents the per-cycle energy dissipated by magnetic hysteresis and eddy currents, better known as core losses. At the design-point stroke of 20 mm (i.e. 10 mm displacement amplitude) the core loss was calculated to be approximately 22 W. In Fig. 7 the open-circuit magnet force is depicted as a surface plot of both displacement and velocity. The force profile can be viewed as evolving in a clockwise ‘path' around the displacement and velocity plane, as is depicted by the arrows. This force map was generated from a bi-cubic interpolation of the FEA simulation results. The surface geometry was then resampled and stored as a three-dimensional array data-set. Then, in the overarching numerical simulation model, the piston's displacement and velocity profiles were cross-referenced to the open-circuit force map using a binary search algorithm with bilinear interpolation, to ultimately reproduce the force-displacement paths shown in Fig. 6. Furthermore, Fig. 7 shows the envisaged 20 mm stroke operating path with the force profile experienced along it. By following the piston motion along the two-dimensional state-space of displacement and velocity, the resulting open-circuit magnetic force experienced by the piston magnet train could more realistically be approximated over using the magnetostatic force profile as a function of displacement only – as is often done in the published literature [59]. It should be noted that this characterising force map is inherently inclusive of the force effects resulting of magnetic hysteresis in the generator core. Consequently, the use of this map to approximate the opencircuit magnetic force constituent allows for the non-linearised magnetic force profile to be accurately emulated, even across changing stroke paths.

4. Finite element analysis of the linear generator The anticipated operation of the linear electric generator was simulated using the transient finite element analysis (FEA) software ANSYS Maxwell 16.0. In this section, simulated results of open and closed-circuited operation is used to infer (map) analytical relations for flux linkage, magnetic force and flux density. As the piston stroke amplitude remains unclear during passive engine operation, the linear generator was characterised across differing piston strokes. This was achieved by subjecting the magnet train to a sinusoidal velocity boundary condition in which the reciprocating frequency was retained at 30 Hz, but the displacement amplitude varied between 2.5 and 15 mm.

4.1. Open-circuited operation The open-circuit coil flux linkage is depicted in Fig. 5 as a function of magnet train displacement. The flux linkage varies approximately linearly with magnet train displacement, despite flattening out at the stroke limits. However, in the engine prototype given in Fig. 1, the power piston's displacement amplitude is in fact mechanically constrained to −12 mm by the engine main body. Any displacement exceeding this limit leads to piston-to-wall collisions, as will be seen from the experimental results in Section 6. The figure further shows that with an increase in magnet train displacement, the flux linkage curves widens. This is due to magnetic hysteresis in the stator cores. Nevertheless, for the range of magnet train displacements simulated, there appears to be a strong correlation between the flux linkage and the magnet train displacement. If the flux widening by magnetic hysteresis is neglected, the average relation between flux linkage and displacement can be reasonably approximated by the fifth-order polynomial:

Fig. 5. Open-circuit flux linkage as a function of magnet train displacement, showing flux widening due to magnetic hysteresis. 9

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

through the generator coil induces an electromagnetic force that, in turn, opposes the movement of the magnet train. When this force is added to the open-circuit magnetic force Fm from Fig. 6, the total electromagnetic force Fmp exerted on the magnet train may be determined. Fig. 10 illustrates the electromagnetic load force divided by load current as a function of magnet train displacement. This load-force profile was determined by subtracting the magnetic force computed from the open-circuit simulation from the total magnetic force computed from the closed-circuit simulation – thus representing only the force constituent induced by current flow. Dividing the electromagnetic load force by the induced current normalises the results obtained from simulating different magnet train displacements, as seen from Fig. 10. The normalised load-force curve can then be approximated by the polynomial: Fig. 6. Open-circuit magnetic force as a function of magnet train displacement, showing force widening due to core loss effects.

F / i (x )= − 1.91 × 10−7x p6 + 9.92 × 10−8x p5 + 1.84 × 10−4x p 4 − 3.14 × 10−5x p3 − 7.33 × 10−2x p2 + 1.59 × 10−3x p + 8.012

(25)

From Fig. 10 it is seen that the simulated data deviates from the fitted polynomial curve at small displacement amplitudes. This deviation is attributed to the magnetic field H not sufficiently penetrating the outer ends of the soft-iron stator during piston reciprocation and thereby inducing only a small magnetic flux path around the outer stator. This results in the straying of the force-per-current curve at the peak displacements, as is clearly evident from the figure. However, nearer to the design-specified stroke of 20 mm, the analytical relation appears to emulate the simulation data sufficiently well. By dividing Eq. (25) by the generator coil length l , the average magnetic flux density B is computed as a function of magnet displacement. Then, from Eq. (20), the Lorentz electromagnetic force is ultimately computed.

Fig. 7. Open-circuit magnetic force as a function of magnet train displacement and velocity, showing the 20 mm path.

4.2. Closed-circuited operation

5. Experimental setup

To investigate the effect of large current flows on generator performance, a purely resistive 1.0 Ω load was simulated as being attached across the coil terminals. The load was simulated by defining an external circuit in Maxwell's Circuit Editor and imposing the circuit as a coil boundary condition. The generator coil winding has an internal resistance of 0.4 Ω that was also accounted for by this boundary condition. Fig. 8 depicts the closed-circuit coil flux linkage as a function of magnet train displacement. In contrast to the open-circuit case that was depicted in Fig. 5, the closed-circuit flux linkage paths appear more spread out. This additional flux widening is attributed to the electrical field around the coil that has to be overcome in order for the flux linkage polarity to be reversed. This increased flux widening induces additional hysteresis losses in the generator core that has not been taken into account in this study. Between the design-point displacement amplitudes of −10 and 10 mm, the average gradient in flux linkage remain approximately linear and similar to that of Eq. (24). For this reason, the assumption of the approximately static coil flux linkage, given by Eq. (24), has been retained. In a follow-up study, for improved accuracy, a characterising map of flux linkage may be derived as has been done for the open-circuit magnetic force. In Fig. 9 the induced load current is depicted for different piston strokes. As seen, the current is asymmetrical about the zero displacement Y-axis. This is primarily due to magnetic hysteresis and internal coil inductance (i.e. generator inductance). This phase difference in the current, that becomes more significant at larger piston strokes, may be reduced by employing a series tuning capacitor in the load circuit, as was shown in Fig. 3 for example. According to Eq. (20), the current flow

To verify and validate of the numerical model, in this section, two experimental test setups are presented. Firstly, only the linear generator is considered. By forcefully reciprocating the power piston assembly, the generator FEA model was validated and the generator characteristic equations verified. Secondly, the fully-assembled engine is considered. By passively operating the assembled engine, the complete numerical model was validated.

Fig. 8. Closed-circuit coil flux linkage as a function of magnet train displacement for different piston strokes, showing flux widening due to the increase in electrical field. 10

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

frequency drive, of which the driving frequency was varied between 10 and 50 Hz – correlating to a maximum reciprocator frequency of approximately 25 Hz on the four-pole motor. As the crankshaft rotates, the connecting rod moves the piston through a fixed 20 mm stroke (10 mm displacement amplitude). The experimentally measured generator terminal voltage and load voltages profiles were used to validate the finite element model and subsequently, indirectly verify the generator characteristic functions derived in Section 4. Experimental tests were performed first with no generator load (open-circuited condition) and thereafter with a purely resistive load consisting of two 2 Ω, 50 W resistors (Arcol HS50) connected in parallel to form a 1 Ω load (more precisely 1.08 Ω±0.6%). In each of these tests the terminal voltage VT , g was recorded using the HBM Quantum MX840A data logger's 60 V analog-to-digital converter (ADC) with sensor accuracy of ±0.5%. In the case of generator loaded operation, the real electrical power delivered to the purely resistive load was calculated. Using an Escort ELC-131D LCR meter, the 70-turn generator coil (consisting of AWG-14 magnet wire) resistance was measured to be 0.4 Ω ± 0.6% and its inductance to be 2.2 mH. By electrically exciting the generator coil with a 10 V, 30 Hz sinusoidal signal, the coil inductance was found to induce a voltage-current phase difference of only 16°. Since this phase difference allows for a good power factor of 96%, a tuning capacitor was omitted from the load circuit. Following a similar approach, the 170-turn motor coil (AWG-19 magnet wire) was measured to have an internal resistance of 1.4 Ω ± 0.6% and inductance of 30.1 mH. Since in this study the displacer motor was used only to kickstart the engine into self-sustained operation, no tuning capacitor was applied in the motor circuit.

Fig. 9. Load current as a function of magnet train displacement.

5.2. Experimental setup of the assembled engine The experimental setup of the fully-assembled engine is depicted in Fig. 12 (left) along with the main equipment used. The engine and its ‘tripod' support stand were bolted to a steel test bench to adsorb and dissipate any engine vibrations. The heating load was supplied by a 1.0 kWth ceramic heating jacket, whereas the cooling water was recirculated to the water reservoir depicted. The assembled engine was purged and pressurised with Grade 5.0 helium via the high pressure line shown. A safety relief valve ensures that the engine cannot be pressurised to above 2.0 MPa gauge. The displacer motor terminals were connected to a sinusoidal voltage generator VT , m to kick-start the engine into operation, whereas the generator terminals were connected to the

Fig. 10. Load force per current as a function of magnet displacement.

5.1. Experimental setup of the linear generator To evaluate the performance of the linear electric generator, the main body with the assembled linear generator was bolted to a stationary support base. As shown in Fig. 11, a three-phase, four-pole induction motor drives the crank. The motor is powered by a variable

Fig. 11. Photograph of the piston linear reciprocator (left) used to validate the generator finite element analysis. Quarter-sectioned rendering of the linear generator as assembled to the reciprocator (right) with the main components as labelled. 11

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Fig. 12. Photograph of the assembled prototype and test measurement instruments (left). Quarter-sectioned rendering of the assembled engine (right), depicting the measurement sensors with the variable names and the locations where they were positioned at.

resistive load with range 1–7 Ω. Using the RCL meter, the electrical wiring to the resistive load bank was measured to have an electrical resistance in the range of 0.16 Ω. In the numerical model, this additional resistance was added to the 0.4 Ω generator coil resistance. An HBM QuantumX MX840A data acquisition computer (DAQ) was used to record the gas pressures, generator voltages and currents at a sampling rate of 2.4 kHz. The Quantum further acted as a half-bridge amplifier to the strain gauge sensors – used to infer the displacer and piston motions from the respective flexure springs. Strain measurements were also recorded at a sampling frequency of 2.4 kHz, with a digital 120 Hz low-pass Butterworth filter being applied during results post-processing. A National Instruments (NI) cDAQ was used to record the engine temperatures measurements. All measurements were logged, processed and time-synchronised using LabVIEW. Fig. 12 (right) depicts the main components and measurement sensors that make up the test setup, along with the locations where these sensors were positioned. The pressure variation in the expansion space Pe was measured using a temperature-compensated, zero to 2.5 MPa (gauge) Wika S-20 pressure transducer with cooling fins, mounted at the top of the engine. A similar pressure transducer was used to measure the backspace pressure variation, Pbs , at the bottom of the engine. Both transducers are accurate to 0.5% of their full-scale readings. However, as the temperature of the working fluid at the thin-film dome is likely to be in the range of 30–80 °C, an additional 1.0% of span measurement error can be expected due to temperature effects. If sensor drift and signal hysteresis are neglected, the root squared error percentage in terms of measurement span is given [60]:

Tot. error =

(0.5)2 + (1.0)2 = 1.18% of span

temperature Tk (between the cooling jacket bottom and the pressure casing flange) and backspace wall temperature Tbs were also measured using T-type thermocouples. In-between the heater and cooler sections, three K-type thermocouples were used to measure the thermal gradient in the engine casing wall, with measurements denoted Tr , top , Tr , Tr , bot . The generator terminal voltage was sampled using the Quantum's 60 V ADC, whereas a 20 A, 75 mV shunt resistor (0.2% accuracy) was used to infer the load current, sampled using the Quantum's 100 mV ADC. To infer the displacer and piston motions, half-bridged strain gauges (HBM type LY11) were attached to opposite ends of the respective flexure springs. By measuring the deflection in the flexure springs, the displacer and piston motions were inferred. The gauges were positioned at the outer legs of the flexure springs where the maximum spring bending deformation occurs. The sensors were calibrated for flexure mid-point displacements ranging between −13 mm and +13 mm on an MTS Criterion 445 testing machine. The Criterion's cross-head displacement with a measurement accuracy of 0.5% was used to calibrate the measured strain reading to displacement. To do so, the strain gauge sensitivity readings were recorded as mV/V measurements with displacement increments of 0.5 mm. Over the anticipated 10 mm stroke, the measurement sensitivity would range between zero and 0.7 mV/V. During this calibration process, a total of 30 sensitivity samples were recorded at each spring displacement point. Results post-processing pointed to a negligible standard deviation and consequently, the sensors were deemed as accurate as the MTS cross-head measurements were (i.e. ±0.5% accurate). In the final engine assembly, all internal electrical wires were passed through three removable pressure-tight feed-through plugs, located at the bottom of engine. This includes the strain gauge sensor wiring with signals x d and x p and the generator and motor terminal wires with signals denoted VT , g and VT , m . All sensor wires were passed underneath the test table before a safety cage was installed over the entire assembly.

(26)

This relates to an error band of ±30 kPa for both measurements. Ktype thermocouples were used to measure the temperature gradient across the heater head, denoted Th, top , Th, mid and Th, bot with an accuracy of ±1.5 °C. The electrical heating jacket with supply voltage G was controlled in a binary-state (on/off) to maintain an average heater head set-point temperature of Th, set . At the cold end, the cooling jacket inlet and outlet temperatures Tw, i and Tw, o were recorded using T-type thermocouples (accuracy of ±0.5 °C). The cooling water was recirculated to a 200 L reservoir at a fixed volume flow rate of 5 L/min. The cold-end

5.3. Experimental procedure The experimental procedure for testing of the fully-assembled prototype may briefly be outlined as follows: 12

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Fig. 14 depicts the open-circuited and closed circuit load voltage across different reciprocation frequencies, whereas in Fig. 15, the peak electrical power delivered to the load is shown. Power trend lines are drawn to visually depict the relationship between voltage and frequency between zero and 50 Hz. Both the open-circuit and closed-circuit results are seen to be in good agreement with the simulated results. For reciprocating frequencies above 20 Hz, there appears to be an increasing over-prediction in the simulated open-circuit voltage results. Due to significant mechanical vibration, the experimental testing was limited to a maximum reciprocation frequency of 25 Hz. Consequently, the performance of the generator could only be estimated at 30 Hz from the trend lines. Fig. 15 depicts how the simulated and measured load power PL increases exponentially with increasing frequency. The FEA model predicts a peak generator output power of 120 We at the 30 Hz design-point frequency, whereas from the experimental trend line, a peak electrical power output of approximately 100 W is estimated. Note that this is nearly four times the 30 We power produced at 15 Hz, which should highlight the importance of increasing the engine operating frequency, should a larger power generator be envisaged for future application. Table 3 summarises the main findings of this experiment, including the deviation percentage between the model and experiment. The good agreement between the simulated and experimental results, even for differing operating frequencies, clearly validates the generator FEA model. This implies that the analytical relations for flux linkage, magnetic force and magnetic flux density derived throughout Section 4 are verified, albeit indirectly.

1. Before engine start-up, assembled engine is pressurised with helium to the specified charge pressure. In this paper, 1.80 MPa (abs.) was used as a demonstrative case study. 2. The heater set-point temperature is specified. The engine is left unoperated until the average temperature profile in the metal heater head reaches the set-point value. In this study, the heater was set to 600 °C. 3. The displacer motor is briefly excited with a 30 Hz, 20 V sinusoidal excitation voltage that ‘kick-starts' the engine into operation. The displacer motor is switched off when the engine becomes operational, and the engine is left to operate passively without any feedback control. 4. The generator load is manually adjusted from 1 Ω to 7 Ω. In this paper, only 1 Ω loaded operation is discussed in detail as a case study. No load or displacer feedback control was investigated in this study. 5. During engine shutdown, the electrical heater is first switched off. The engine is left operating in a generator loaded condition until stalling occurs. Cooling water is recirculated until the average hotside temperature decreases to below 100 °C. The engine was operated under normal atmospheric conditions near sea-level in a lab controlled environment. The water coolant flow rate was maintained at 5 L/min with an inlet temperature of approximately 31 °C. 6. Results and discussion Prior to operating the assembled engine, the main engine components were analysed to identify and verify principal modelling assumptions and engine characteristic parameters. This includes experimentally verifying the spring stiffnesses, verifying the displacer and piston first mode natural frequencies and calculating the mechanical friction forces of the assembled engine. The numerical model was then validated in two different experimental setups. Firstly, to validate the linear generator FEA model, the power piston magnet train was mechanically reciprocated as per the test setup shown in Fig. 11. Secondly, the engine was fully assembled as per the experimental setup depicted in Fig. 12. The engine was test-operated to hot-end temperatures ranging from 400 to 600 °C and charge pressures of between 1.0 and 1.80 MPa absolute. However, in this study, only a selected experimental case study of realistic engine operation is discussed. Experimental results are qualitatively compared to results obtained from the simulation model, to verify and validate the numerical model, and to investigate and quantify model deviations.

6.2. Validation of self-sustained engine operation Self-sustained engine operation becomes possible when the differential gas pressure forces acting on the displacer shaft are sufficiently large to overcome the frictional and gas shuttling losses over the engine heat exchanger channel [61]. As a case study of realistic and envisaged engine operation, in this experiment, the engine was pressurised to 1.80 MPa absolute with helium and operated to a heater set-point temperature of 600 °C. Following the experimental procedure outlined in Section 5.3, self-sustained operation was achieved for nearly six minutes at a time before engine stalling occurred due to a decline in hot-end temperature. To verify that the numerical model is capable of emulating the operational behavioural characteristics of the FPSE prototype, a qualitative comparison between simulation and experimental results is made for an illustrative case study of self-sustained operation. Table 2 lists the operating parameters and engine temperatures recorded for a selected case study of 1 Ω generator loaded operation. Direct comparisons are

6.1. Validation of the generator FEA model In this experiment, the linear reciprocator shown in Fig. 11 was used to forcefully reciprocate the power piston assembly over a fixed 20 mm stroke. Tests were performed with reciprocation frequencies varying from 5 to 25 Hz. Subsequently, experimental measurements of both open- and closed-circuited operation were compared to numerical results obtained from the linear generator FEA model. Fig. 13 compares the open-circuit voltage VOC and the resultant closed-circuit load voltage VL of the 1.08 Ω load being attached across the generator terminals, when the piston was reciprocated at the peak permissible frequency of 25 Hz. The modelled results are seen to be in good agreement with the experiment, despite the open-circuit voltage amplitude being slightly over-predicted by 2 V from the experiment. This deviation is due to unaccounted-for core loss effects that play a prominent role in opposing the magnetic flux density in the open-circuit operation. Unlike the simulated results, the experimental voltages appear slightly more asymmetrical, suggesting that the magnet train was in fact not reciprocated perfectly about the centre of the stator core during testing, likely due to a small offset in the experimental setup.

Fig. 13. Comparison between simulated and measured open- and closed-circuit generator voltage and current. 13

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Fig. 14. Simulated and experimentally measured open- and closed-circuit voltages at different reciprocation frequencies.

Fig. 16. Simulated and measured engine metal temperature profile.

6.2.1. Dynamic steady-state operation Fig. 16 depicts the simulated casing wall Twll , bore cylinder wall Tcyl and regenerator mesh temperatures Tr , mesh as a functions of distance x from the heater top. The figure gives both the initialised (Tinit ) temperature profile at simulation start and the final model-converged temperatures at simulated dynamic steady-state operation. The figure serves to compare the experimentally measured casing wall profile Twll, exp . to the simulated results. As seen, the heater head wall temperature decreases nearly linearly towards the regenerator side, arguably due to large axial conductive heat transfer in the engine metal walls. In fact, with reference to the energy balance diagram in Fig. 25, the simulation results suggest that approximately 230 W of heat flows from the heater wall to the cooler wall via the engine pressure casing, and that additionally, 200 W is lost to cooling jacket via the displacer's bore cylinder wall. Overall, a good correlation between the measured and simulated temperature profile is noted. In Fig. 16 the modelled regenerator mesh temperatures have the appearance of a ‘temperature band' on the scale of this figure, with the individual mesh cells fluctuating by approximately +5 °C along an otherwise linear temperature profile across the regenerator section. At the cold end, where the water cooling jacket is located, the casing wall temperature is maintained below 100 °C. The model predicts a 2.61 °C increase in cooling water temperature, which is in good agreement to the experimentally determined 2.9 °C increase. Due to conductive heat transfer occurring across the width of the internal heating and cooling fins, the bore cylinder wall temperature Tcyl is similar to that of the engine casing wall. However, this thermal contact acts to increase the axial conductive heat losses from heater to cooler, leading to the additional 200 W of heat being rejected to the water. Fig. 17 depicts the displacer and piston displacements relative to their mid-stroke locations. The model predicts self-sustained operation at an engine frequency of 29.97 Hz, with peak-to-peak displacer and piston strokes of 26.2 mm and 25.3 mm respectively, and a pistondisplacer phase difference of 69°. This is in good agreement to the measured displacer and piston strokes of 25.8 mm and 25.9 mm respectively. However, the experiment shows the power piston to collide with the engine's main body cylinder once during every second compression stroke at a displacement amplitude of −12 mm. In contrast, the model does not predict collisions to occur due to being slightly under-predictive of the piston stroke. As a result of these metal-onmetal collisions ‘bouncing' the piston away from the main body wall, the experimental piston-displacer phase lag varies from 68° just before an impeding collision, to 55° just after a collision. The slight deviation between measured and simulated piston phasing should be noted from Fig. 17 as a result thereof. In comparing the engine displacements with the pressure variation given in Fig. 18 – as illustrated by the locations of

Fig. 15. Simulated and experimentally measured peak electrical power delivered to 1.08 Ω load at different reciprocation frequencies. Table 2 Measured and simulated engine operating parameters for a demonstrative case study of self-sustained engine operation. Test case

Model input

Unit

No. of cells per heat exchanger Set-point temperature Abs. charge pressure, P Generator load, Rload Heater head top, Th, top

– 600 1.80 1.18 576.5

73 – 1.80 1.18 576.5

– °C MPa Ω °C

Heater head bottom, Th, bot Regen. top-wall, Tr , top

521.2 473.3

521.2 –

°C °C

Regen. mid-wall, Tr Regen. bot-wall, Tr , bot Cold-end wall, Tk Backspace wall, Tbs Coolant inlet, Tw, i Coolant outlet, Tw, o Coolant flow rate, fw

307.1 117.1 61.3 34.2 30.7 33.6 5.0

– – – – 30 – 5.0

°C °C °C °C °C °C 1/min

made between the simulated and measured casing wall temperature profile, the displacer and piston displacements, the expansion space and backspace pressure swings, and the generator output voltage and currents. The simulated flow of heat through the engine is analysed to investigate the engine's theoretical thermal performance. An energy balance of simulated engine operation concludes the results discussion.

14

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Fig. 19. Expansion and compression space pressure versus chamber volume.

Fig. 17. Displacer and piston displacement over time, showing the occurrence of a piston collision and its influence on the piston's phase.

To reaffirm whether the model is over-predicting the pressure variation – as has been assumed in this study – or whether the experimental pressure swing is in fact slightly under-measured, would require additional testing using more sensitive pressure transducers. Fig. 19 depicts the per-cycle variation in the expansion and compression space pressures relative to their respective chamber volumes, with the arrows depicting how the pressure curves evolve with volume. In this pressure volume (PV) diagram, the gas pressure in the compression space has been assumed to be equal to that measured in the expansion space (i.e. the gas pressure drop in the regenerator section has remained unmeasured in this experiment). By integrating the area encapsulated by the simulated PV curves it is calculated that, during a single engine cycle, the gas has 9.58 J of work done onto it in the compression space, whilst the gas delivers 16.08 J of positive work in the expansion space. This relates to boundary work rates of 287.09 W and 481.83 W respectively, with the workspace delivering a net positive indicated power of 194.53 W. Not shown here is that the model predicts 0.16 W of net work being done onto the backspace gas by the reciprocating piston and displacer, which is virtually negligible. From the experiment is calculated that, over a two cycle average, the engine delivers 343.47 W of indicated power in expanding the working fluid, whereas 178.42 W of indicated power is required to compress the working fluid. This leads to a net positive indicated power of 165.05 W. The difference in experimental and simulated indicated power is largely attributed to the over-prediction in workspace pressure by the numerical model, but also due to the piston collisions that influence the displacer-piston phase difference, thereby leading to the cyclical pressure swings and lower indicated powers. Fig. 20 shows the expansion and compression space pressures as functions of the total workspace volume. The simulated pressures in the respective chambers are nearly indistinguishable from each other, apart from viscous dissipation losses in the regenerator channel, and a small phase shift between the expansion and compression space pressures that results due to the finite sonic speed at which the pressure between these spaces are balanced out. The encapsulating area within the simulated expansion space curve indicates a net work rate of 212.61 W, whereas the compression space curve yields 214.42 W. The resulting net difference of 1.82 W between the indicated powers is attributed to a slight difference in phase that exist between the gas pressures within the chambers. When compared to the simulated net indicated power of Fig. 19, it is noted that the average net indicated power is 18.99 W larger as what was previously calculated. This difference arises as the indicated power calculated using the total engine volume diagram does not consider the work being dissipated to overcome internal flow friction. This finding also highlights a limitation of the experimental test setup: since only the expansion space pressure was measured, the

Fig. 18. Expansion space and backspace pressure over time, showing the decrease in pressure swing as induced by the piston collision.

maximum and minimum pressure swings – it is calculated that peak simulated pressure swing leads the displacer motion by 48° and the piston by 117°. This is in good agreement with the average experimental phase differences of 50° and 110°. The effect of the piston-casing collision is also clearly noticeable from the engine gas pressure variation shown in Fig. 18, for which the expansion space pressure swing are seen to be increasing-and-decreasing periodically. This periodic variation in pressure swing is a direct result of the periodic variation in displacer-piston phase difference, further also correlating to a periodical variation in the extracted boundary work. In addition to the collision-induced mismatch in simulated and measured pressure swing, the numerical model over-predicts the pressure swing peak by nearly 80 kPa. Due to the complexity of the overall numerical model, this over-prediction may be attributed to a range of factors including (but not being limited to): lower-thansimulated heat transfer occurring in the heater section, unaccounted-for deadspace in the expansion space manifold and in the pressure transducer cooling channel, larger-than-simulated fluid dissipation losses, and the simplifying assumptions of adiabatic expansion and compression. As a counter argument to the pressure variation being over-predicted by the model, it may also be argued that the experimental pressure variation is under-measured instead. Arguably, this could be as a result of the expansion space transducer exhibiting an insufficient response settling time to accurately characterising the pressure peaks.

15

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

̇ ̇ . The cooler jacket to the engine casing wall, denoted Qheat and Qcool figure also shows the net heat transfer rates to and from the working fluid, Qḣ − g and Q̇g − k . These heat transfer rates are compared to the simulated electrical power being delivered to the generator load, Pgen . The model predicts that, in order for the engine to be operated in thermodynamic steady-state at the hot-end temperature profile specified in Fig. 16, the electrical heater is required to deliver 1212.3 W of thermal power. However, since the existing electrical heater was limited to 1000 W, it is postulated that the decrease in heater temperature that led to eventual engine stalling after six minutes of operation is attributed to this approximately 200 W deficit in heating input power. At the cold end, 998.6 W of heat is rejected to the cooling water and additionally, 19.3 W of heat is rejected by the backspace casing (not shown here). The net rate of heat input to the engine may be calculated to be 194.4 W, which is also equal to the per-cycle average mechanical power delivered to the piston and displacer. From Fig. 23 the cyclical variation in net heat absorption rate by the hot gas and net heat rejection rate by the gas in the cooler section should also be noted. The two distinctive peaks in the heater-to-gas and gas-to-cooler heat exchange rates occur when the bulk of the working fluid is being shuttled through the heat exchanger channel. When steady-state operation is achieved, the net difference of the time-integrated area under the heater and cooler rates ̇ ̇ ) is equal to the net indicated power being delivered by (Qheat and Qcool the working fluid minus the heat rejected by the backspace. In Fig. 24, the net heat transferred to the hot gas, Qh − g , and the net heat rejected by the cold gas, Qg − k , is given to highlight the scale of heat transfer occurring within each of the respective heat exchanger sections. When thermodynamic steady-state is reached, the net amount of heat transferred from the working fluid to the regenerator mesh metal during a single engine cycle is zero. The sinusoidal appearance of Qg − rm illustrates the thermal charging and discharging of the regenerator mesh . It must be noted that the regenerator mesh stores and delivers nearly five times the heat energy delivered by the hot-end. As the regenerator surface area and volume are increased, this ratio will increase further and so will the thermodynamic efficiency, but at the expense of larger regenerator frictional pressure losses and dead volume. Consequently, the development of an optimised regenerator presents itself as a trade-off study between increasing the regenerator thermal capacity (and dead volume) and minimising frictional pressure drop [62]. The significantly large amount of unused heat being rejected to the cold-end should also noted. The net difference between Qh − g and Qg − k represents the work converted to mechanical energy minus the heat rejected by the backspace.

Fig. 20. Expansion and compression space pressure versus total engine volume.

calculated experimental indicated power of 165.05 W remains inclusive of viscous dissipation losses in the regenerator section. An additional pressure sensor, located at the compression space side, is recommended for future study should it be required that the indicated power be more accurately measured or the work lost to viscous dissipation in the regenerator channel be studied. Fig. 21 compares the simulated emf, load voltage and current to the experimentally measured load voltage and current. The simulated load voltage amplitude is nearly 7 V lower than that of the induced generator emf e . Due to the inductive nature of the generator coil, the simulated and measured load voltages lag behind the induced emf by 11.87°. Overall, the simulated load voltage and generator current appear to be in good agreement with that of the experiment, despite the model over-predicting the load voltage and current amplitudes by 2.4 V and 2.1 A respectively. These over-predictions of the numerical model are mainly attributed to the simplifying open-circuit relation used to model the generator flux linkage, i.e. Eq. (24), that does not take into account the widening or flattening of the flux linkage curve under closed-circuited generator operation. Around 60.07 seconds the experimental current deviates in phase from the simulated results, correlating to just after the piston had collided with the engine body. The experimentally measured and simulated electrical output powers are compared in Fig. 22. It is seen that model generally overpredicts the peak in load power by between 30 and 40 W. This is as a result of the over-predictions in load voltage and current previously described. The cyclic variation that is seen in the simulated load power profile is a result of the piston being offset slightly towards the positive displacement side (into the backspace). In contrast, the experimental load power clearly indicates a linear ‘ramp-up’ in output power, from approximately 60 W up to 70 W over a period of four power strokes. This ramping-up period indicates that the piston collides with the engine casing only during every second engine cycle, correlating to the periodic decrease in workspace pressure swing in Fig. 18 that is seen to occur every second engine cycle. Arguably, this ramping up of the output power outweighs the otherwise anticipated cyclic variation in power that is shown by the model. When averaged over the duration of two piston stroke cycles, the experimental engine is calculated to deliver an average electrical power of 24.1 W, which is twenty percent lower than the model's predicted aveerage power of 31.4 W. But in spite of the model's over-prediction in electrical power output, the overall shapes of the modelled and measured power profiles are in reasonable agreement.

6.2.3. Energy balance of simulated operation An energy balance diagram of simulated engine operation is given

6.2.2. Simulated heat rate and energy transferred To aid in analysing the engine's theoretical thermal performance, Fig. 23 gives the net simulated heat transfer rates from heater and

Fig. 21. Generator load voltage (V) and current (I) over time. 16

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

50% the heat rejected to the cooling jacket. The figure indicates that 235.33 W of heat is conducted axially from the heater towards the regenerator in the engine pressure casing wall alone. Since there is nearly any heat transfer between the regenerator wall and working gas (only 2.13 W is convected to the gaseous working fluid), the bulk of this axial parasitic heat loss is rejected by the cooling jacket. If the internally located heating and cooling fins are both in contact with the casing and bore-cylinder walls, a thermal short-circuit is induced. This further results in 264.74 W of heat being conducted from the heater, across the fin breadth and into the displacer bore cylinder, of which 247.74 W is rejected to the cooling water. Fig. 25 also illustrates how the workspace and backspace gas reserves are interlinked with each other through gas leakage and pistondisplacer boundary work. By gas enthalpy leakage, energy is transported at a net rate of 19.14 W from the compression space towards the backspace. The working space exerts boundary work at a rate of 192.56 W onto the power piston and 1.97 W onto the displacer. In turn, the free-running power piston exerts 1.28 W of boundary work onto backspace, of which 1.12 W is exerted back onto the displacer shaft. The displacer dissipates mechanical energy at a net rate of 3.09 W to contact-sliding mechanical friction. Due to the large piston mechanical friction, a considerable amount of 114.07 W of mechanical power is expended to piston-sliding friction. As a result of the poor electricalconversion efficiency of the linear generator, 45.79 W is irreversibly dissipated by generator core and ohmic losses, so that only 31.42 W is ultimately converted to usable electrical power. Finally, from the energy balance in Fig. 25, the net thermal efficiency may be calculated to be 16.03%; the mechanical efficiency to be 39.72%; and the electrical efficiency to be 40.69%. This leads to a considerably low net thermal-to-electrical efficiency of only 2.59%. It is suggested that the low net efficiency is attributed to the significant axial conductive heat losses in the casing and bore cylinder walls that, in effect, limit the available heat input to expanding the workspace gas to 712.23W. The importance of reducing this large parasitic loss in a future study should be emphasised. Additionally, large mechanical friction forces dissipate nearly sixty percent of the indicated power being produced. Should a more idealised case study be presented, one in which the axial heat losses are neglected, the thermal efficiency could theoretically be increased to 192.56/(712.23)=27.03%. Furthermore, should the combined mechanical and electrical efficiencies be assumed

Fig. 22. Power delivered to the electrical load over time.

Fig. 23. Simulated heat input and heat rejection rates over two cycles, compared to the generator electrical power in Watts.

Fig. 24. Simulated heat transferred to and from the respective heat exchanger channels over two cycles, given Joules.

in Fig. 25, to visually illustrate the flow of energy through the engine. It should be noted that a significant percentage of the thermal heat input is rejected by the cooling jacket. In fact, 82% of the heater input power, or nearly five times the energy converted to mechanical power, is rejected to the cooling water. This poor thermal performance is mainly attributed to the significant axial conductive heat transfer that occurs through the casing and bore cylinder walls, making up approximately

Fig. 25. Energy balance diagram of simulated engine operation, depicting the per-cycle averaged rates of energy transport in Watts. 17

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

as 90%, which would entail a more optimised laminated generator core with reduced ohmic losses (and possibly incorporating dynamic gas bearings to reduce piston friction and improve its alignment), the still unoptimised prototype could be expected to deliver a net efficiency of nearly 24%. To further increase the net efficiency, a detailed parametric optimisation study on engine design parameters could be conducted using the validated numerical model.

•

6.3. Summary of model and test results Table 3 summarises the measured and simulated outcomes of both experiments. The deviation percentage, defined as the percentage deviation of the numerical results to the experiment, is given in an effort to quantify the model's agreement with the experiment. It should be noted that, although the piston reciprocation experiment was performed at a lower frequency and stroke, both the model and test results show a slightly larger electrical power being produced than that of selfsustained operation. This decrease in load power is attributed to the additional 0.16 Ω electrical resistance of the generator load bank that was not present in the piston reciprocation test. The table shows that the modelled pressure swing deviates by almost twenty percent from the experiment. In addition to the piston-metal collisions that reduced the experimental indicated power, the over-estimation in modelled pressure swing contributes to a deviation of 15.15% in workspace indicated power. Furthermore, despite the over-predictions in pressure and indicated power, the simulated displacer and piston motions deviate by less than five percent, the piston-displacer phase difference by less than 13%, and the pressure-piston phase by less than 7%. The prototype engine delivered on a two-cycle average 24.10 W of electrical power, whereas the model predicted 31.42 W. Consequently, the overall model over-predicts the engine output power by 26.30%. Still, this should be considered a remarkable agreement when considering how complex and comprehensive the overall numerical model really is.

•

For future work, three follow-up studies can be recommended:

• Firstly, pertaining to the existing engine prototype, it is deemed

critical to reduce the large axial heat losses in the engine metal walls. This would in theory allow the engine thermal efficiency to be increased beyond sixteen percent. The improvement of the mechanical and electrical efficiencies are also recommended, but

7. Conclusions and recommendations This paper has given an overview of the work done in modelling, manufacturing and experimentally testing a unique free-piston type Stirling engine (FPSE) electric generator. In presenting a uniquely comprehensive and validated free-piston Stirling engine model that, for the first time, computes the engine thermodynamic and kinematic behaviour in an explicitly coupled manner, this study has made a significant contribution to the still-limited body of knowledge surrounding the modelling of free-piston Stirling engine generators. The main conclusions are drawn briefly as follows:

Table 3 Summary of simulated and experimentally measured results with the percentage deviation given in square brackets.

• The paper firstly gave an overview of the design and manufacture of

•

characteristic equations, realistic linear generator performance could be simulated without explicitly invoking a coupled finite element analysis. Thirdly, the paper presented an experimental investigation in which measured results were qualitively compared to their numerically simulated equivalents, in order to validate the numerical model. The linear generator model was first validated before a demonstrative case study of engine operation was presented and discussed in detail. The paper found that the numerical model over-predicted the piston by 2.61%; the piston-displacer phase difference was overpredicted by 12.42%; the workspace indicated power was overpredicted by 15.15%; and the average electrical power output was over-predicted by 23.30%. Moreover, the deviations in piston phasing and indicated power were mainly attributed to physical piston-casing collisions that occurred during engine operation, although an additional investigation into pressure over-prediction by the model was recommended. The study concludes that the coupled thermodynamic-kinematic model discussed in this paper is indeed capable of emulating the dynamic operational behaviour of a real FPSE from first principles. As no such comprehensive and validated numerical model of a FPSE is currently available from the published literature, this study presents a significant contribution to the still-limited body of literature surrounding the modelling of free-piston Stirling engine operational mechanics. It is envisaged that this validated numerical model will greatly aid prospective engine designers with the modelling, analysis and ground-up design of more improved free-piston Stirling engines in general.

a 100 We (peak) FPSE prototype. The geometric layout of the engine prototype has been based on Sunpower Inc's P2A engine and Microgen's commercially available 1 kW engine, but with the unique addition of a linear motor that is directly attached to the displacer shaft. A schematic overview of this novel, laboratory-based engine prototype was presented and the nominal operating and design parameters were given. The experimental setup was outlined with the measurement sensors and experimental variables clearly defined. Secondly, the paper gave an overview of a third-order numerical model entitled SENSE (Stirling Engine Numerical Simulation Environment) in which the transport equations of the engine's compressible gaseous working fluid remain inherently coupled to the engine kinematics and electrodynamics. A transient electromagnetic finite element analysis of the linear generator was performed so that simplifying analytical relations for flux linkage, magnetic force and magnetic flux density could be mapped against magnet displacement and velocity. In the unique modelling approach outlined by this study, by employing these mapped 18

Test

Model

Unit

[% Dev]

Piston reciprocation Piston stroke Reciprocation frequency Open-circuit voltage (p-p) Closed-circuit voltage (p-p) Load current (p-p) Peak load power (p-p) Average load power Mechanical power required Electrical efficiency

20.00 25.00 22.25 17.19 15.92 68.42 26.16 – –

20.00 25.0 25.87 18.29 16.94 77.46 31.68 59.18 53.53

mm Hz mm V A W W W %

[–] [–] [13.99] [6.01] [6.02] [11.67] [17.42] [–] [–]

Self-sustained operation Displacer stroke Piston stroke Operating frequency Piston-displacer phase diff. Pressure-piston phase diff. Pressure-displacer phase diff. Pressure swing (p-p) Load voltage (p-p) Load current (p-p) Heat input rate Cooler heat rejection rate Exp. space indicated power Comp. space indicated power Workspace indicated power Avg. electrical power out Net engine efficiency

25.83 25.95 29.98 60.21 110.08 49.87 472.24 18.20 15.24 – 1012.34 343.47 178.42 165.05 24.10 <2.4

26.23 25.29 29.97 68.75 117.32 48.57 577.63 22.54 18.94 1212.25 998.61 481.83 287.09 194.53 31.42 2.59

mm mm Hz ° ° ° Pa V A W W W W W W %

[1.52] [−2.61] [−0.03] [12.42] [6.17] [−2.68] [18.25] [19.25] [19.54] [–] [−1.37] [28.72] [37.85] [15.15] [23.30] [7.34]

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

•

•

would constitute to a complete redesign of the linear generator. In a revised design, the linear generator could be made to induce a nonlinear increasing magnetic force so that piston over-stroking be limited. To decrease piston mechanical friction, the magnetic sideload forces acting between the magnet train and the generator core should be minimised. The alignment of especially the power piston could be improved by employing a dynamic gas bearing. Secondly, it is strongly recommended that an active feedback controller be developed so as to ensure the continual operation of the engine without collisions occurring. Different active control strategies should be investigated using the numerical model and laboratory-based prototype – this includes the investigation of direct displacer control, or active piston load control. For future work it is recommended that the as-developed numerical model first be used to test and investigate different control strategies directly. Thereafter, with the insight gained from the modelling results, an engine controller should be designed, implemented and tested. Thirdly, it is recommended to conduct a comprehensive parametric study on engine geometry so that a more optimised and task-suited engine can ultimately be developed. It is recommended that the design of a larger output engine should start with the development of a more powerful and optimised linear generator. Only when the linear generator has been designed and characterised should the numerical simulation model, presented in this paper, be employed as a design tool to further investigate and refine the engine geometry.

Latham, New York, U.S.A.; 1999. [13] Goldwater B, Piller S, Rauch J, Cella A. Demonstration of a free piston stirling engine driven linear alternator, Tech. Rep. EY-76-C-2764. Latham, New York: Mechanical Technology Incorporated; 1977. [14] Schreiber JG, Geng SM. RE 1000 free-piston Stirling engine hydraulic output system description. Cleveland, Ohio U.S.A: National Aeronautics and Space Administration; 1987. [15] Oriti S, Scott W. Advanced Stirling convertor (ASC-E2) performance testing at NASA glenn research cent, National Aeronautics and Space Administration (NASA) Glenn Research Center, Cleveland, Ohio, U.S.A.; 2011. [16] Wilson SD, Schifer NA, Williams ZD, Metscher JF. Overview of Stirling technology research at NASA Glenn Research Center, Technical Publication NASA/TM—2016218909. U.S.A.: National Aeronautics and Space Administration (NASA), Glenn Research Center, Ohio; 2016. [17] Gibson M, Oleson S, Poston D, McClure P. NASA's kilopower reactor development and the path to higher power missions, techreport NASA/TM-2017-219467. U.S.A.: National Aeronautics and Space Administration (NASA) Glenn Research Center, Ohio; 2017. [18] Qiu S, Gao Y, Rinker G, Yanaga K. Development of an advanced free-piston stirling engine for micro combined heating and power application. Appl Energy 2019;235:987–1000. [19] CSP World, The end for dish-Stirling CSP Technology? Infinia files for bankruptcy, [Online]. Available at: http://www.csp-world.com/news/20130929/001196/enddish-stirling-csp-technology-infinia-files-bankruptcy, [2015, March 3]; 2013. [20] Microgen, Microgen Engine Corporation, [Online]. Available at: http://www.microgen-engine.com/, [2018, April 1]; 2018. [21] Qnergy, Smart, efficient, reliable energy solutions with the stirling engine, [Online]. Available at: https://www.qnergy.com/, [2018, February 6]; 2018. [22] ARPA-E, sustainable economic MCHP stirling (SEMS) generator, [Online]. Available at: https://arpa-e.energy.gov/?q=slick-sheet-project/sustainable-economic-mchpstirling-sems-generator, [2019, April 4]; 2019. [23] Zare S, Tavakolpour-Saleh A. Free piston Stirling engines: A review. Int J Energy Res 2019:1–32. [24] Ramos JA. Thermodynamic analysis of stirling engine systems, Ph.D. thesis, KTH Royal Institute of Technology, School of Industrial Engineering and Management, Sweden; 2015. [25] Ahmadi MH, Ahmadi M-A, Pourfayaz F. Thermal models for analysis of performance of stirling engine: a review. Renew Sustain Energy Rev 2017;68:168–84. [26] Walker G, Senft J. Lecture notes in engineering: free piston stirling engines. 1st ed. Heidelberg, Germany: Springer-Verslag; 1985. [27] Araoz JA, Cardozo E, Salomon M, Alejo L, Fransson TH. Development and validation of a thermodynamic model for the performance analysis of a gamma stirling engine prototype. Appl Therm Eng 2015;83:16–30. [28] Alfarawi S, Al-Dadah R, Mahmoud S. Enhanced thermodynamic modelling of a gamma-type stirling engine. Appl Therm Eng 2016;106:1380–90. [29] Martini WR. Stirling engine design manual, national aeronautics and space administration (NASA), Cleveland, Ohio, U.S.A., 2nd ed.; 1983. [30] Snyman H. Second order analysis methods for stirling engine design, master's thesis, Faculty of Engineering, University of Stellenbosch, South Africa; 2007. [31] Redlich RW, Berchowitz D. Linear dynamics of free-piston stirling engines. Proceedings of the institution of mechanical engineers, Vol. 199, Ohio, U.S.A. 1985. p. 203–13. [32] Karabulut H. Dynamic analysis of a free piston stirling engine working with closed and open thermodynamic cycles. Renew Energy 2011;36:1704–9. [33] Formosa F. Coupled thermodynamic-dynamic semi-analytical model of free piston stirling engines. Energy Convers Manage 2011;52(5):2098–109. [34] Strauss JM. Direct piston displacement control of free-piston stirling engines, Ph.D. thesis, Faculty of Engineering, University of Stellenbosch, South Africa; 2013. [35] Tavakolpour-Saleh A, Zare S, Bahreman H. A novel active free piston stirling engine: modelling, development, and experiment. Appl Energy 2017;199:400–15. [36] Dyson R, Wilson S, Tew R. Review of computational stirling analysis methods. 2nd international energy conversion engineering conference. 2004. [37] Urieli I. A computer simulation of stirling cycle machines, Ph.D. thesis, Faculty of Engieering, University of the Witwatersrand, South Africa; 1977. [38] Urieli I, Berchowitz D. Stirling cycle engine analysis (modern energy studies). 1st ed. CRC Press; 1984. [39] Schock A. Stirling engine nodal analysis program. J Energy 1978;2(6):354–62. [40] Goldberg LF. A state space and continuum mechanics analysis of stirling cycle machines, Ph.D. thesis, Faculty of Engineering, University of the Witwatersrand, Johannesburg, South Africa; 1987. [41] Goldberg LF. One- and two-dimensional stirling machine simulation using experimentally generated reversing flow turubulence models, Tech. Rep. NASA-CR185285, National Aeronautics and Space Administration (NASA) Lewis Research Center, Ohio, U.S.A.; 1990. [42] Yuan S, Spradley I. A third order computer model for stirling refrigerators. Adv Cryog Eng 1992;37:1055–62. [43] Andersen SK. Numerical simulation of cyclic thermodynamic processes, Ph.D. thesis, Technical University of Denmark, Department of Mechanical Engineering, Denmark; 2006. [44] Deetlefs IN. Design, simulation, manufacture and testing of a free-piston stirling engine, Master's thesis, Faculty of Engineering, University of Stellenbosch, South Africa; 2014. [45] Wang K, Dubey S, Choo FH, Duan F. A transient one-dimensional numerical model for kinetic stirling engine. Appl Energy 2016;183:775–90. [46] Shunmin Z, Guoyao Y, Jongmin O, Tao X, Zhanghua W, Wei D, et al. Modeling and experimental investigation of a free-piston stirling engine-based micro-combined

CRediT authorship contribution statement B.J.G. de la Bat: Conceptualization, Data Curation, Validation, Writing - original draft. R.T. Dobson: Conceptualization., Writing review & editing T.M. Harms: Writing - review & editing. A.J. Bell: Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] Hudson H. From rural village to global village: telecommunications for development in the information age. 1st ed. New York, U.S.A: Taylor and Francis; 2013. [2] Murugan S, Horáak B. A review of micro combined heat and power systems for residential applications. Renew Sustain Energy Rev 2016;64:144–62. [3] Skorek-Osikowska A, Remiorz L, Bartela Lukasz, Kotowicz J. Potential for the use of micro-cogeneration prosumer systems based on the stirling engine with an example in the polish market. Energy 2017;133:46–61. [4] Organ AJ. Thermodynamics and gas dynamics of the Stirling cycle machine. Cambridgeshire, U.K.: Cambridge University Press; 1992. [5] Beale WT. Free piston stirling engines - some model tests and simulations, SAE International. [6] Senda FM. Aspects of waste heat recovery and utilisation (WHRU) in pebble bed modular reactor (PBMR) technology, master's thesis, faculty of engineering. South Africa: University of Stellenbosch; 2012. [7] Renzi M, Brandoni C. Study and application of a regenerative stirling cogeneration device based on biomass combustion. Appl Therm Eng 2014;67(1):341–51. [8] Dochat GR. Design study of a 15 kW free-piston stirling engine - linear alternator for dispersed solar electric power systems, techreport DOE/NASA/0056-79/1. Latham, New York, U.S.A: Mechanical Technology Incorporated; 1979. [9] Caballero GEC, Mendoza LS, Martinez AM, Silva EE, Melian VR, Venturini OJ, et al. Optimization of a dish Stirling system working with DIR-type receiver using multiobjective techniques. Appl Energy 2017;204:271–86. [10] Gibson M, Mason L, Bowman C, Poston D, McClure P, Creasy J, et al. Development of NASA’s small fission power system for science and human exploration, techreport NASA/TM-2015-218460. Ohio U.S.A.: National Aeronautics and Space Administration (NASA) Glenn Research Center; 2015. [11] Gibson MA, Poston DI, McClure P, Godfroy TJ, Briggs MH, Sanzi JL. Kilopower Reactor Using Stirling TechnologY (KRUSTY) nuclear ground test results and lessons learned, technical report NASA/TM-2018-219941. Ohio U.S.A.: National Aeronautics and Space Administration (NASA) Glenn Research Center; 2018. [12] Dhar M. Stirling space engine program, final report, Mechanical Technology Inc.,

19

Applied Energy 263 (2020) 114585

B.J.G. de la Bat, et al.

Lewis Research Center, Ohio; 1996. [55] Çengel YA, Ghajar AJ. Heat and mass transfer: fundamentals and applications, 4th ed. New York, U.S.A: McGraw-Hill; 2011. [56] Kröger D. Air-cooled heat exchangers and cooling towers. 1st ed. Stellenbosch, South Africa: University of Stellenbosch; 1998. [57] Boldea I, Nasar SA. Linear electric actuators and generators. Cambridge, U.K.: Cambridge University Press; 1997. [58] Lenand DK, Priest JF, Keiter DE, Schreiber JG. Development of a power electronics controller for the advanced stirling radioisotope generator, technical publication NASA/TM-2008-215031. U.S.A.: National Aeronautics and Space Administration (NASA) Glenn Research Center, Cleveland, Ohio; 2008. [59] Kim J-M, Choi J-Y, Lee K-S, Lee S-H. Design and analysis of linear oscillatory singlephase permanent magnet generator for free-piston stirling engine systems. AIP Adv 2017;7(5). 056667. [60] Hashemian HM, Jiang Jin. Pressure transmitter accuracy, ISA Transactions 48 (4) (2009) 383–8. [61] Berchowitz D. Operational characteristics of free-piston stirling engines. In: 23rd intersociety energy conversion engineering conference, Vol. 1, Denver, Colorado, U. S.A.; 1988. p. 107–12. [62] Ibrahim MB, Tew RC. Stirling converter regenerators. Florida, U.S.A.: CRC Press: Taylor & Francis Group; 2012.

heat and power system. Appl Energy 2018;226:522–33. [47] Gedeon D. Sage user's manual, Gedeon Associates, Athens, Ohio, U.S.A, 11th Edition; 2016. [48] Geng S, Mason L, Dyson R, Penswick L. Overview of multikilowatt free-piston stirling power conversion research at Glenn Research Center, Tech. Rep. NASA/ TM—2008-215061. Ohio, U.S.A.: National Aeronautics and Space Administration (NASA) Glenn Research Center; 2008. [49] Peckham engineering and tool, precision linear flexure bearing cartridge, Philips Laboratory, Kirtland Air Force Base, Albuquerque, New Mexico, U.S.A.; 1994. [50] B. J. G. de la Bat, Theoretical simulation, manufacture and experimental evaluation of a free-piston stirling engine electric generator, Ph.D. thesis, Faculty of Engineering, University of Stellenbosch; 2019. [51] Patankar S, Spalding D. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int J Heat Mass Transf 1972;15(10):1787–806. [52] Bird R, Stewart W, Lightfoot E. Transport phenomena. 2nd ed. New York, U.S.A: John Wiley and Sons Inc.; 2002. [53] Grinnell S. Flow of a compressible fluid in a thin passage, transactions of ASME: Paper No. 55-SA-13; 1956. p. 765–76. [54] Gedeon D, Wood J. Oscillating-flow regenerator test rig: hardware and theory with derived correlations for screen and felts, technical publication NASA TECHREPORT-198442. U.S.A.: National Aeronautics and Space Administration (NASA)

20