Control-Oriented Modeling of the Dynamics of Stirling Engine Regenerators

Stirling devices are energy conversion devices, which can be operated as heat engines or as heat pumps. In engine mode, they can operate using any heat source such as external combustion, waste heat, or solar thermal power. Although the basic Stirling engine design has been around for over two centuries, newer varieties continue to be periodically invented, with the free piston Stirling engine [1] being the most widely known modern example. Due to their ability to use any heat source, there has been a resurgent interest in using Stirling engines for utilization of solar thermal power [2–4], as well as in micro-combined heat and power applications [5–8] (with at least one recent commercial offering by Whispergen^®, which produces household scale, Stirling engine-based micro-combined heat and power units). On the research side, there has also been a recent surge of interest in the dynamics and control of Stirling engines [5,9–16]. Given the important role that electronic controls have had in improving the performance of the internal combustion engine in the past few decades [17], it is arguable that control engineering can have a significant impact on performance of analogously modern versions of the Stirling engine.

Of particular interest to the authors is the new concept of the actively controlled Stirling engine [9,12] in which actuator control of the displacer piston allows for significant freedom in the design of the engine's limit cycle. Potential improvements in both power density and efficiency can be substantial [9].

Control-oriented dynamic models based on the Schmidt assumptions have been shown, via experimental validation, to be able to accurately predict the operating frequency of a free piston Stirling engine [11,14]. However, higher fidelity models may be necessary for the purposes of engine limit cycle design in an actively controlled engine where the displacer is actuated by a control system. For example, the Schmidt assumptions do not capture the dynamics of the gas states, which neglects the finite heat transfer rates within the engine. If used for limit cycle design, the Schmidt assumptions could result in a suggested operating frequency that is faster than the true optimal, by not allowing sufficient time for heat to be exchanged between the gas and the engine materials. To this end, we present a control-oriented model, which captures the dynamics of the working gas states, placing particular emphasis on modeling arguably the most important component of a Stirling engine, namely, the regenerator. As a historical aside, Stirling's original patent was not for the air engine (which was known at the time), but rather for the regenerator [18] (which he termed, the “economizer”).

The regenerator acts as a type of thermal capacitor that stores and retrieves thermal energy to the working gas during different phases of the engine's cycle. A detailed description of its operation is provided in Sec. 2. It is made up of a porous matrix of high thermal capacitance material through which the working gas flows. The working gas and the matrix interact via convective heat transfer. Thus, a model of these dynamics involves the equations of compressible flow together with convective heat transfer through porous media. The partial differential equations (PDE) modeling these processes are computationally intensive to simulate and difficult to use for control design. We therefore develop a hierarchy of reduced models in this paper, which ultimately yields low-order models with high fidelity for the typical frequencies of cyclic operation.

There is a significant challenge that arises in the present model reduction problem, which is also relevant in a wider context. The regenerator operates in feedback with the dynamics of the remaining parts of the engine, which are assumed to be lumped. Thus, the reduced model needs to approximate not only the regenerator's open-loop behavior, but more importantly, the behavior of the overall coupled system. This is an instance of closed-loop (or “feedback-aware”) model reduction, which arises not only in the problem of controller reduction [19–22] but also in other contexts such as large-scale model reduction with the requirement to preserve certain properties such as passivity [23–25]. In this paper, we argue that since the regenerator will ultimately operate periodically when coupled to the remainder of the engine, a chirp-proper orthogonal decomposition (chirp-POD) technique (a POD technique with a chirp signal as the input) over the frequency range of potential operation provides an effective method for this nonlinear model reduction problem.

The paper is organized as follows: In Sec. 2, we first introduce and summarize the basic operation of the regenerator in a Stirling engine. We then develop a one-dimensional PDE model using the Euler equations of gas dynamics coupled with gas/mesh heat transfer. There are four fields in this coupled system of PDEs, three for the gas state (our choice is density, velocity, and pressure) and an additional field for the mesh matrix temperature. In Sec. 3, these equations are nondimensionalized to uncover four nondimensional parameters that indicate potential time-scale separations. While one of those parameters is the well-known Mach number, the other three are lesser known and are specific to convective heat transfer between an oscillating gas and a mesh material. These parameters characterize various types of possible engines, and we pay special attention to this characterization. In Sec. 4, we develop a hierarchy of singularly perturbed models by taking the respective small parameter limits, as well as averaging the mesh dynamics under the additional assumption of periodic operation. It turns out that the most useful and accurate model for typical operation is the one where the fast acoustic dynamics are removed, and we reduce that particular model further using a chirp-POD technique in Sec. 5. Section 6 concludes the paper about the utility of these model reduction techniques in active control design for Stirling engines.

Figure 1 is a diagram of a generic Stirling engine. A sealed working gas shuttles back and forth through the regenerator between the cold and hot sections, each section being in thermal contact with external cold and hot heat exchangers, respectively. This movement of the gas is primarily caused by the displacer piston, and while the power piston's motion has some effect on the gas movement, it is a smaller effect than that due to the displacer, and can be neglected for a preliminary understanding of the engine's operation. The role of the regenerator will be explained below, but it can temporarily be thought of as simply a low-pressure-loss connection between the hot and cold sections.

A Stirling engine is a special type of air engine. The basic operating principle of an air engine is simple. At any one time, the hot and cold sections of the engine are at roughly the same pressure, and that pressure is dependent on the average temperature across those sections. As the working gas shuttles back and forth between the two sections, more or less of its volume is contained in either the hot or cold sections, and therefore the average temperature oscillates, which, in turn, causes the pressure to oscillate. These pressure oscillations then drive the working piston to perform mechanical work on a load.

The cycle described previously is driven by the displacer's movement. In traditional Stirling engines, displacer motion is induced through kinematic linkages from the power piston. Little power is required to drive the displacer since the working gases at either side of it are roughly at the same pressure. Linkages designed so that displacer and power pistons are approximately $90 deg$ out of phase will typically produce a stable oscillation [26]. Alternatively, in free piston Stirling engines [1], the displacer and power piston are dynamically (rather than kinematically) linked through a gas spring (the so-called bounce space). Figure 1 shows these kinematic/dynamic linkages only conceptually since they differ from one type of Stirling engine to another.

A more recent concept is that of the actively controlled Stirling engine [9,12], where the displacer's motion is directly actuated by a control actuator whose motion can be designed to optimize the engine [9].

One of the motivations of the present work is that a more realistic active control design will probably require a better dynamical model than the commonly used idealized Schmidt model. An additional goal is to develop a modeling framework for a generic Stirling engine, i.e., largely independent of the presence or type of piston/displacer linkages. We therefore emphasize the most dynamically complex, and arguably the most important part of the engine, that is the regenerator.

The regenerator is an open connection between the cold and hot sections, which is filled with a metal (or graphite) mesh of material. It acts like a porous channel for gas flow, but yet significant convective heat exchange occurs between the flowing gas and the mesh. The purpose of the regenerator is to act as a “thermal capacitor.” For example, without the regenerator, as hot gas flows into the cold section, most of its excess heat will be rejected to the cold exchanger. The regenerator retains some of that thermal energy in the mesh material, and the gas enters the cold section colder than it would have otherwise been without a regenerator. In the other half of the cycle, the gas flowing back to the hot section retrieves some of the heat stored in the mesh, and arrives at a higher temperature than it would otherwise be at (without the regenerator), thus needing to extract less thermal energy from the hot exchanger to reach the hot side temperature. This concept of a thermal capacitor was the basis of Stirling's original patent [18], which is referred to as the “Economizer.” Without the regenerator, an air engine would have very low efficiency and power output.

We begin this section with a distributed one-dimensional PDE model of the regenerator based on the Euler equations of gas dynamics together with mesh/gas heat exchange. We then describe a standard lumped model for the thermodynamics of each gas section, and finally the dynamics and possible kinematic linkages of the pistons are described. These three separate pieces of the model are schematically shown in Fig. 2.

The regenerator consists of channel through which the working gas flows back and forth between the hot and cold gas section. This channel is typically filled with a mesh (a metal or graphite) material, which runs lengthwise down the tube. The gas and the mesh material exchange heat as the working gas flows back and forth through the channel. The main feature of our model of the regenerator is that due to gas flow being mainly in the axial direction, all variables are well approximated as being constant along axes perpendicular to the flow. For the gas, this leads to the Euler equations of one-dimensional gas dynamics, while the mesh material's temperature can be modeled by the one-dimensional heat equation. The thermal interaction between the mesh and the gas is captured by a simple model of convective heat transfer.

The geometry of the one-dimensional model is depicted in Fig. 3. The dynamics of the gas are given by the one-dimensional Euler equations, which reflect the conservation of mass, momentum, and energy, respectively,

where ρ, v, p, E, T, and $Φ$ are the spatially distributed gas' density, velocity, pressure, energy, temperature, and mesh material temperature, respectively. All of these fields are functions of space x and time t, which are suppressed for notational simplicity. The terms in the last column represent the nonconservative effects of viscous friction and mesh–gas heat exchange, respectively. The latter term represents a simple Fourier law of heat exchange between the gas and mesh with a heat transfer coefficient of k_g. It is a simplification which combines the effects of conductive and convective heat transfer in a single velocity-independent coefficient [13].

where c_m, ρ_m, and k_m are the (specific) heat capacity, density, and heat conduction coefficient of the mesh material, respectively. This coupling between the gas dynamics and the mesh thermal state is depicted in Fig. 4.

Thus, there are four possible combinations of boundary conditions depending on whether the gas is flowing in or out of each end of the regenerator. This “state-dependent switching” of boundary conditions presents a special challenge in modeling and numerical approximations of this system.

The method we use to simulate such systems is the essentially nonoscillatory (ENO) scheme, which is commonly used to for compressible flow equations. The ENO scheme will not be covered here, but interested readers are referred to Refs. [27–29]. An important feature of the ENO scheme is that it is able to account for the state-dependent switching of the boundary conditions. It does so by using an “upwinding” scheme, where the direction of upwinding switches with the sign of the velocity. Thus, the ENO scheme is able to capture this state-dependent switching automatically.

The first is the ideal gas law, and the second is the expression for the gas internal thermal energy.

where a₀ and a₁ are the cross-sectional areas of the boundaries between the hot and cold sections, respectively, and the regenerator. $v(.)$ and $ρ(.)$ are the velocity and density fields as determined by the gas dynamics model of the regenerator. Note the switching nature of the above equations. If the velocity v(0) at the boundary between the hot section and the regenerator is positive, then the hot section is losing mass of density ρ₀, the density of the mass already in the hot section. If v(0) is negative, then it gains mass of density $ρ(0)$⁠, i.e., the density of the gas at the left boundary of the regenerator.

which is equal to the rate at which (thermal + kinetic) energy is flowing in/out of the section, plus the work rate of pressure forces, minus the mechanical work done by the gas on its surroundings due to volume changes, and the last term represents the conductive heat transfer (with coefficient K_i) between the heat exchanger at the walls of the section, which are assumed kept at constant temperatures H_i. Note that T(i), $ρ(i)$⁠, and v(i) are the boundary values of fields determined by the distributed regenerator model described in Sec. 2.1.1.

We finally note the number of state variables and inputs. Both volume and volume change rate V_i and $V˙i$ can be regarded as inputs imposed by the pistons' positions and velocities. The boundary velocities v(i) can also be regarded as inputs imposed by the regenerator model. The remaining variables T_i, ρ_i, P_i, and E_i are dynamic state variables that are constrained by two (for each section) static algebraic relations (9), and therefore one can choose only two state variables for each lumped section. The goal was to express the fact the both of the two pistons and the regenerator's all have their own dynamics. This is summarized in Fig. 2.

where v(i), $ρ(i)$⁠, V_i, and $V˙i$ are regarded as inputs.

where x_p and m_p are the position and mass of the power piston, respectively, $P1−Pex$ is the pressure difference between the cold section and the external side of the power piston, $cpx˙p$ is a viscous dissipation term, which models the power dissipated in the load, x_d is the displacer's position whose velocity is assumed directly assignable by an external control input u. The time-varying section volumes V₀ and V₁ and their derivatives $V˙0$ and $V˙1$ can be considered as outputs of this dynamical system as shown above (the equations for $V˙0 and V˙1$ are not shown) with a_d and a_p the pistons' cross-sectional areas, and $V¯0 and V¯1$ as the respective volumes when x_p = 0, x_d = 0.

In beta-type kinematically linked Stirling engines, there are further algebraic constraints. The displacer and power piston are connected to a flywheel, which provides the feedback necessary for a stable limit cycle to form. The kinematics for the beta engine can be expressed using the geometrical relations from Fig. 5 as follows:

where I and θ are the moment of inertia and angular position of the flywheel, respectively, F_p is the reaction force between the power piston and the flywheel, $ϕ$ is the phase difference between the two pistons, r_p and r_d are the radial attachment locations of the pistons on the flywheel. The equations above assume that the displacer and arms connecting the pistons to the flywheel are massless. The latter are also assumed to be sufficiently long so that the forces they exert on the flywheel and pistons are essentially horizontal.

Note that although there are four states in Eqs. (14) and (15), the two algebraic constraints (16) and (17) reduce the number of state variables in this model to two. We choose θ and $θ˙$ as coordinates for these states.

The parameters used in the sequel for numerical examples for both the mechanical and the section models are based on the work of Mehdizadeh and Stouffs [10], where a martini type Stirling engine is modeled. The engine modeled here is a simplified version of that engine. The choice of wall temperatures, helium as the working gas, approximate dimensions, and nominal pressure are all taken from this source.

where $ρ¯ and p¯$ are nominal gas density and pressure, $Φ¯$ is nominal mesh temperature, ω_b and $v¯$ are the engine's frequency (T_b is the engine's period) and nominal advection speed, respectively. The latter two parameters are typically not known in advance, but can be estimated based on other parameters and the engine's geometry. Note that in the nondimensional variables, the ideal gas law becomes simply $p̃=ρ̃T̃$⁠.

ε₁ is the ratio of the advection time (average time it takes the gas to advect from one end of the regenerator to the other) to the engine's period. This ratio is less than one for a well-designed engine; otherwise, motion reversal in the regenerator occurs before gas has shuttled from one section to the other. The Mach number ε₂ is the smallest of the previous parameters and is typically much less than one. A reasonably efficient engine would also have a much smaller head loss (through the regenerator) than the nominal operating pressure, and therefore ε_f is typically much less than one. Finally, the mesh material is designed to have a high thermal capacity, and therefore the mesh thermal time constant $cpρm/kg$ is typically much longer than the engine's period T_b.

It is now clear that the first three equations are candidates for a singular perturbation type of model reduction while the last equation may be simplified using an averaging technique. These simplifications are detailed in Sec. 4. We close this section with a further discussion on which properties of a particular engine determine how small each of the parameters is.

where V_r is the void volume in the regenerator, and V_s is the volume displaced by the displacer during one stroke of the engine.

Given these approximations, it is possible to justify the assumption that all of the epsilons are small. If the flow speed in the regenerator does not become supersonic, which should be the case for most Stirling engines, then $ε2=(v¯/c¯)$ will be less than one. If the void volume (V_r) is less than the stroke volume $(Vs)$⁠, then $ε1=(Vr/πkVs)$ must be less than one. This should also be the case for most Stirling engines as a well-designed Stirling engine minimizes the dead volume (volume that is not part of the expansion or compression process) throughout the engine, which the regenerator volume is considered to be. Because of the large thermal inertia of the matrix material, the engine frequency is much faster than the matrix temperature dynamics. Therefore, $(Tb/(cpρm/kg))$ will be small. The final small parameter $εf=(βLv¯/p¯)$ should also be less than one since $βLv¯=Δp$ where $Δp$ is pressure difference across the regenerator. Most Stirling engine assumptions, including the common Schmidt assumptions, assume the pressure difference across the regenerator is negligible compared to the nominal pressure. Thus, it is safe to assume that for most Stirling engines, these values will all be less than or much less than one.

In this section, we develop a hierarchy of singularly perturbed models in progressing order of coarseness. The first model (42) is concerned with only the gas' states, but eliminates acoustic dynamics. This is applicable to low Mach number conditions, which is typical for Stirling engines. The second model involves (Bogoliubov) averaging of the mesh matrix dynamics motivated by its relatively large thermal inertia and periodic operation. This leads to the conclusion that the solid mesh temperature profile is linear (50). We also present a quasi-steady state (QSS) model that parametrizes the temperature, pressure, and density distributions across the regenerator with a small number of parameters that are algebraically related to the boundary conditions.

The above equations represent flow, compression, and thermal exchange with the mesh at the advective time-scale. The much faster acoustic dynamics have been removed from this model.

which follows from $2(ppx)x=(p2)xx$⁠. We note that the dynamics of p are primarily diffusive with the temperature difference $(Φ−T)$ as a distributed source term. The dynamics of ρ on the other hand are advective (with the advection velocity $−(1/εf)px=v$⁠, and a regenerative term proportional to p_xx).

The pressure equation in Eq. (42) is of the diffusive type with boundary conditions given at both ends. It can therefore be easily discretized using a central difference scheme. The density equation, however, is of the advection type with switching boundary conditions according to Eq. (8) (note that $v(x,t)=−px(x,t)/εf$⁠, and therefore the switching boundary conditions on ρ depend on p in the no-acoustics model). We therefore use a standard first-order upwinding scheme to discretize ρ_x, where at each grid point, the direction of upwinding is based on the local velocity $−px(x,t)/εf$⁠.

where $Φ0$ and $Φ1$ are constants of integration determined by the boundary conditions. In our case, the ends of the mesh are assumed to be in thermal contact with the section walls and as such, the boundary conditions are that the temperature at the ends of the mesh must be equal to the temperature of the section walls. As was mentioned at the start of this section, averaging analysis indicates that the difference between this and the true limit cycle is of order ε_m. Therefore, the matrix material will be assumed to be a fixed linear profile, which interpolates between the two wall temperatures.

As an aside, we note that the linear temperature profile (50) would have also been a consequence of assuming that the thermal conductivity of the mesh (k_m) is much greater than the combined conduction and convection coefficient between the gas and the mesh (k_g). However, a large $km/kg$ essentially implies a “thermal short” between the hot and cold exchangers through the regenerator, which leads to low efficiency. Therefore, a well-designed engine would not have a large $km/kg$ ratio, and we did not invoke this assumption here.

In engines where the product $ε1εf$ can be considered small (e.g., low head loss through regenerator and/or slower engine frequency than advection time through regenerator), a significant simplification occurs in the equations. It turns out that the regenerator states depend algebraically on the boundary conditions. We term this the QSS model since it represents a situation where the regenerator dynamics are much faster than the remaining engine's dynamics.

where the second equality follows from Eq. (38).

where we have also listed the pressure boundary conditions as well.

The details of this derivation are in Appendix  B.

for some function f, which is a combination of polynomials and exponentials in α. Therefore, given the values of ρ₀, p₀, and p₁, a root finding routine can be used at each time step to solve for the corresponding α.

while the temperature follows from the ideal gas law.

A comparison was carried out of the above QSS model and the full model of Sec. 2 using the ENO scheme. For each of the models, the gas sections had two states each and the flywheel dynamics had two states giving a total of six states for the lumped portion of the engine. The ENO scheme for the full model used 100 grid points while the QSS model of the regenerator is memoryless. Therefore, the total state count is $4·100+6=406$ for the full model and 6 for the QSS model, giving a significant reduction in state dimension.

For comparison, both beta engine models were given the same initial condition and simulated until a steady-state limit cycle was reached. Figure 6 shows a comparison of the limit cycles of the sections' states (the profiles of the fields in the regenerator are not shown). The trajectories show that the QSS model is very close to the full model for small values $ε1εf$ while still being qualitatively close for moderate values of $ε1εf$⁠.

The results showed that the simplified QSS model performed similarly to that of the finite difference model as long as the product $ε1εf$ was small. A natural question would be if a higher-order perturbation approximation would yield a better match. This is not the case for this application. It turns out that higher-order perturbation approximations require that the time derivatives of the inputs to the perturbed model (boundary conditions) be known. The calculation of the time derivatives of the inputs (the section states) requires the use of the lower-order approximations of the regenerator profiles. As such, these time derivatives will differ from those of the unreduced system. Via experimentation, it was discovered that the higher-order profile approximations are very sensitive to errors in these time derivatives. As a result, the higher-order approximation performed worse than the lower-order approximation. In the case that the inputs to this regenerator model are not states but instead predetermined time histories whose time derivatives are known, then a higher-order approximation of these profiles would likely yield better results. However, this was not the case for this application.

In this section, we investigate a POD numerical method for model reduction of the no-acoustics model of Sec. 4.1. The main idea behind our particular POD method is in the choice of simulation conditions used to obtain the time traces from which POD modes are extracted. Since the regenerator will operate in a time-periodic manner once connected to other components of the engine, our approach is to generate snapshots from a simulation where all signals follow a chirp profile (a sinusoid with a linearly time-varying frequency) with a frequency range representative of the engine's potential operating frequency range. This approach has some commonality with that used in Refs. [31,32], and we term it a “chirp-POD” technique.

where ${tl}$ and ${xk}$ are the time instances and the grid point locations, respectively. The fields $ρ̃$ and $p̃$ are the fluctuations of density and pressure from the nominal solution corresponding to zero displacer motion.

The above procedure is applied to the density and pressure data separately.

where ${ϕkρ}$ and ${ψkp}$ are the first N columns of the corresponding matrix $Φ$ in Eq. (58) expressed as functions of a continuous variable (by, e.g., interpolation). The reduced regenerator model then has 2N states. We will refer to this reduced model as the chirp-POD model.

We now present results of the chirp-POD model reduction technique performed using snapshots from a simulation of the no-acoustics regenerator model (42) connected to the beta-type engine model. In order to generate trajectories with a time-varying instantaneous frequency (i.e., a chirp), an input is needed. This was done in the beta engine model (14)–(19) by imposing a time-varying trajectory ${θ(t)}$ on the flywheel. To do this, the dynamics of the flywheel/power piston assembly (Eqs. (14) and (15)) were removed from the model, and the imposed θ trajectory was used to generate piston motions from Eqs. (16) and (17).

The chirp input consisted of increasing the flywheel's speed linearly in time from 40 to 600 rad/s. This frequency range was chosen to include 80 rad/s, which is the nominal operating frequency of this engine model. It was found that a choice of truncation order of N = 3 produced very good results, while N < 3 resulted in unrealistic trajectories. This assessment was carried out as follows. The chirp-POD modes were obtained from simulations with a forced flywheel trajectory, while the comparison of trajectories was done with dynamics of the beta engine (14)–(19) connected to the regenerator models. These models were simulated until stable limit cycles emerged, and the comparison is displayed in Fig. 7 (left) showing a very close match.

The reduced model was further tested by connecting it with a sections' model that has double the wall heat transfer coefficient compared to that used in obtaining the chirp-POD modes. The results are shown in Fig. 7 (right), and indicate that the technique of using chirped inputs appears effective in capturing POD modes that work for a large frequency range rather than a single operating condition.

In Sec. 5.1, we compared the reduced and full model by comparing trajectories in two operating regimes. The driven Stirling engine model, however, offers the possibility of performing a more comprehensive comparison since it has a natural input, and therefore one can compare input–output behavior. Since the anticipated operation is periodic, it is natural to attempt to quantify this input–output behavior using a notion of frequency response. For nonlinear systems, a full frequency response would require analyzing the full harmonic content input–output periodic pairs at all possible amplitudes. We instead use a simpler notion of a frequency response assuming inputs to be pure sinusoids. The details are described below.

First, we describe the POD step. In the previous subsection, the kinematically linked engine was forced to have a flywheel trajectory, which followed a chirp signal. For the driven engine, the displacer motion is a direct input, and we therefore used the displacer position as a chirped input with the same range as previously stated, 40–600 rad/s. We used the driven model (13) where the displacer velocity was chosen so that the displacer motion is the required chirp signal. A POD analysis was done on the resulting data and a reduced model using the first four modes was chosen. These POD modes are shown in Fig. 8.

The “frequency response” comparison between the full and reduced model was performed as follows: Both models were simulated with the same pure sinusoid as input. The outputs were then periodic signals with several harmonic components. The amplitude of each of those harmonics was found, and this analysis was repeated for a range of input signal frequencies. Figure 9 (top) shows the results where only the first four harmonics of the outputs are shown (higher harmonics' amplitudes were too small to be relevant). As the results show, the first two or three harmonics match up quite well and deviation only begins once the spectral content has dropped by at least an order of magnitude. Because both models are nonlinear, changing the amplitude of the input signal may result in a different frequency response. To test this, the displacer amplitude was increased by 50% and the experiment was repeated with the same POD modes used previously. Those results are shown in Fig. 9 (bottom) and indicate similarly good reduced model fidelity.

It is worth mentioning that this reduced chirp-POD model of the driven engine has a total state dimension of 15 (eight states for the regenerator model, four for the two gas sections, and three for the piston dynamics). This model gives trajectories that are essentially indistinguishable from the full order 406 states model. In addition to a smaller state dimension, the reduced model also runs much faster. For a particular displacer motion, one cycle of the actuated engine was simulated with the original regenerator model, the acoustic removed version, and the reduced chirp-POD model. The simulation time for the original model was 1203 s, the acoustics removed version took 72 s, and the reduced chirp-POD model took 12 s (100 times faster). Most of this time savings comes from the fact that the time stepping routine can now take larger steps, as there are no longer fast acoustic phenomena that require a very small time step to accurately capture.

We have shown that effective and significant model order reduction is possible for Stirling engine regenerators. Depending on the operating regime of the engine, reductions from model orders in the hundreds to single digits are possible. More generally, our main contribution is the development of a methodology for model reduction of systems with compressible gas dynamics. By identifying several nondimensional parameters that characterize different types of engines, we used singular perturbation methods to obtain a hierarchy of simplified models.

The regenerator has distributed dynamics while the remaining sections of the engine have lumped models. We therefore had to address the issue of how one does POD, or balanced truncation type model reduction, with the objective that the reduced model works well when connected with the remainder of the engine. This might be thought of as a feedback-aware model reduction objective. For the Stirling engine, it appears that the chirp-POD technique we used is particularly suited since the overall system operates eventually in a limit cycling mode. We anticipate that this technique might be useful for other cyclically operating devices common in energy harvesting and thermoacoustics. In particular, when the limit cycle itself is to be designed such as in optimal periodic control, reduced models such as the ones we presented can significantly reduce the computational complexity of optimal control calculations while retaining the requisite fidelity in dynamics. This is the intended future application for this work.

Finally, we note that although this paper is primarily concerned with Stirling engine regenerators, much of our model reduction work is likely applicable to the regenerators (stacks) of thermoacoustic engines and heat pumps [33–36]. The underlying physical mechanisms are quite similar, though the operating regimes in parameter space might be different.

This research is partially supported by NSF Grant No. CMMI-1363386. The authors would like to thank Professor Frederic Gibou of the Mechanical Engineering department at UCSB for guidance on numerical methods, especially the ENO scheme.

where we have used Eq. (1) to go from the second line to the third. The rest follows from the chain rule.