Stirling Engine Simulation

Introduction

When I was studying for my exit exams, I couldn’t find a good Stirling engine simulation online. There were plenty of simplified diagrams showing the four stages of the cycle and their corresponding positions on a p–V diagram, but I couldn’t find one that showed the engine and its thermodynamic cycle moving continuously.

So I decided to make one myself, partly because it helped me understand how the engine works, and partly in the hope that it might help future students understand it a little better.

The engine shown here is best described as a beta-type Stirling engine. The design is based on a small low-temperature tabletop Stirling engine I got for Christmas a few years ago. A displacer moves the working gas between a hot and a cold space, while a separate work piston changes the total volume. The two move 90° out of phase.

For the thermodynamics I use the ideal Schmidt model.

Model and Assumptions

The simulation uses the ideal Schmidt model of a Stirling engine. This model simplifies the thermodynamics enough that the pressure can be calculated directly from the positions of the displacer and work piston.

For this simulation I assume that:

the working gas is an ideal gas

the amount of working gas is constant

pressure is uniform throughout the engine at any instant

the hot space remains at a constant temperature T~h

the cold space remains at a constant temperature T~c

piston and displacer motion is sinusoidal

the engine runs at constant angular velocity

pressure losses and mechanical losses are ignored

These assumptions reduce the engine to two gas spaces: a hot space and a cold space. Their volumes change with the crank position, but each space is assigned a fixed temperature and both are assumed to have the same instantaneous pressure.

This also means that gas transferred between the two spaces is effectively assumed to reach the temperature of its new space immediately. R

Simulation

The simulation below is interactive. Change any of the inputs to see how they affect the engine and its p–V cycle.

Cold-side temperature (T꜀) [K]

Hot-side temperature (Tₕ) [K]

Amount of gas (n) [mol]

Rotational frequency (f) [Hz]

Displacer swept volume (V_s,d) [m³]

Work-piston swept volume (V_s,w) [m³]

Hot-side dead volume (Vᴅ,h) [m³]

Cold-side dead volume (Vᴅ,c) [m³]

FPS GIF

ㅤ

Physics

I described the state of the engine in terms of two angles: the crank angle $\theta$ (shown in green in the simulation) and the work piston angle $\alpha$ .

The crank angle is defined as $\theta = \omega t$ , where $\omega$ is the angular velocity and $t$ is time. Since the work piston lags behind the crank by 90°, its angle is given by $\alpha = \theta - 90°$ .

All other physical quantities describing the engine, including the positions of its components, pressure, and volume, can then be expressed as functions of these two angles.

Geometry

I model the motion of both the work piston (purple) and the displacer piston (green) as sinusoidal, rather than reproducing the exact geometry of the connecting rods. Since the crank radius is much smaller than the connecting rod length, the resulting deviation from the sinusoidal approximation is negligible.

The volume displaced by the work piston is given by

V_w(\alpha) = \frac{V_{s,w}}{2} (1-\sin{\alpha}).

The hot-side gas volume is

V_h(\theta) = V_{d,h} + \frac{V_{s,d}}{2}(1-\sin\theta)

while the cold-side gas volume is

V_c(\theta, \alpha) = V_{d,c} + V_w(\alpha) + \frac{V_{s,d}}{2}(1+\sin\theta).

Here, $V_{s,w}$ and $V_{s,d}$ are the swept volumes of the work piston and displacer piston, respectively, while $V_{d,h}$ and $V_{d,c}$ are the hot- and cold-side dead volumes.

Adding the hot- and cold-side volumes gives the total gas volume:

\begin{aligned} V(\alpha) &= V_h(\theta) + V_c(\theta,\alpha)\\ &= V_{s,d} + V_{d,c} + V_{d,h} + V_w(\alpha). \end{aligned}

As you can see the displacer piston does not change the total gas volume; it merely redistributes the gas between the hot and cold regions. Since all other terms are constant, the total gas volume depends only on the position of the work piston, described by $\alpha$ .

Pressure

The pressure calculation is the main thermodynamic part of the simulation.

For an ideal gas, $pV=nRT.$

That equation cannot quite be applied to the entire engine using a single temperature, because part of the gas is assumed to be at $T_h$ and another part at $T_c$ .

The model instead assumes that the pressure is uniform throughout the engine. The amount of gas in each space can therefore be written as

n_h=\frac{pV_h}{RT_h}

and

n_c=\frac{pV_c}{RT_c}.

The total amount of gas is fixed:

n=n_h+n_c.

Substituting the previous expressions gives

n = \frac{p}{R} (\dfrac{V_h}{T_h} + \dfrac{V_c}{T_c})

Rearranging for pressure gives the equation used in the simulation:

p = \frac{nR}{\dfrac{V_h}{T_h} + \dfrac{V_c}{T_c}}

As the displacer moves gas towards the hot end, the effective temperature of the working gas increases and so does its pressure. Moving it towards the cold end does the opposite.

The 90° phase difference between this pressure variation and the motion of the work piston is what produces the p–V loop.

p–V Diagram and Work

The p–V diagram is probably the most useful part of the simulation.

The code loops through one complete cycle, calculates the volume and pressure every 2 degrees, and plots the resulting data.

The work performed by the gas is

dW=p,dV.

Over one complete cycle, work is

W=\oint p\,dV.

This is visualised in the p–V diagram as the area enclosed by the loop. Since the simulation only calculates pressure and volume at discrete crank angles, this integral is evaluated numerically from 1000 calculated points rather than analytically.

Power

Once the work per cycle is known, the average indicated power is straightforward:

\bar P=Wf

where $f$ is the rotational frequency.

The instantaneous power follows from

P = \frac{dW}{dt} = p\frac{dV}{dt}.

Only the work piston changes the total volume. Its volume is

V_w = \frac{V_{s,w}}{2}(1-\sin\alpha),

with

\alpha = \omega t - 90^\circ.

Differentiating with respect to time gives

\frac{dV}{dt} = -\frac{V_{s,w}}{2}\cos(\alpha)\omega.

The instantaneous power is therefore

P= -p\frac{V_{s,w}}{2}\cos(\alpha)\omega

Other Calculations

The simulation also displays a few other quantities, including the mean pressure, compression ratio, swept volume, and dead volume. These remaining calculations are mostly straightforward or performed numerically. You can find their implementations, along with the rest of the simulation’s source code, on GitHub.

Limitations

As listed above the simulation relies on quite a few assumptions most notably that cold and hot spaces remain at constant temperature. This is obviously not the case in the real world as heat transfer takes time and the working gas can not reach the temperature of the hot or cold side instantly.

This becomes particularly important as the engine speed increases. In the model, frequency does not affect the thermodynamic state of the gas: the same pressure and volume are calculated at a given crank angle regardless of how quickly the engine is rotating. The p–V loop, and therefore the work per cycle, consequently remains unchanged with frequency.

Since

\bar P = Wf,

the predicted power therefore increases linearly with frequency. A real Stirling engine does not behave this way. At higher speeds, finite heat-transfer rates, flow restrictions, pressure losses, friction, leakage, and other effects become increasingly important.

The model also does not attempt to describe the detailed behaviour of a regenerator or other non-ideal processes within the engine.

For these reasons, the simulation is not an an accurate predictor of the performance of a real Stirling engine. Instead, it should be used to help understand how these engines work and to accomplish my original goal of visualising the different stages of a Stirling engine cycle. A more realistic simulation would require a considerably more complicated thermodynamic model.

Implementation

I wrote the entire simulation in JavaScript and used the HTML Canvas API for both the animation and the p–V diagram. This means that it can run entirely in the frontend. The calculations required by the thermodynamic model are computationally inexpensive; drawing, updating, and exporting the animation are considerably more demanding.

Animation

The animation is implemented using requestAnimationFrame(), which calls the drawing function before each browser repaint.

Rather than incrementing the crank angle by a fixed amount every frame, I calculate it from the elapsed time and angular velocity:

const theta = (omega * t) % (2 * Math.PI);
const alpha = (theta - Math.PI / 2 + 2 * Math.PI) % (2 * Math.PI);

This allows the animation to run at approximately the correct speed regardless of the browser’s frame rate.

For each frame, the canvas is cleared and the entire engine is redrawn.

The p–V diagram is generated separately by calculating the pressure and volume at 180 equally spaced positions over a complete revolution. These values are only recalculated whenever an engine parameter changes.

This isn’t how I would structure a serious thermodynamic simulation program, but it has the significant advantage that it works.

GIF Export

I also added a button to export the animation as a GIF, mostly because I wanted an easy way to use it elsewhere, for example in a PowerPoint presentation.

For this, I used the JavaScript library modern-gif. The export function renders the animation at a user-defined frame rate and collects the resulting canvas frames. Instead of recording the animation in real time, it calculates the engine’s position separately for each frame.

The collected frames are then passed to modern-gif, which encodes them into a GIF. The resulting file is made available through a temporary browser URL.

One disadvantage is that GIF encoding can be rather slow, particularly because the canvas has a (hard-coded) resolution of 1000 × 1500 pixels. At the default rotational frequency of 0.15 Hz, a complete revolution takes approximately 6.7 seconds, requiring about 100 frames at 15 FPS.

Code

The complete source code is available on GitHub.

There are certainly cleaner ways of implementing quite a lot of it. In particular, the drawing code contains a somewhat excessive number of hard-coded coordinates, and several calculations could probably be consolidated into fewer functions.

I could refactor all of this, but the simulation does what I originally wanted it to do, and I have other things I’d rather spend my time on.