Internal Combustion Engine Simulation in 2D: Pressure, Temperature, and Flow

using Rasmah
using Statistics

Overview

The four-stroke engine turns the slider-crank motion of a piston into a pressure rise that does work. In this tutorial we simulate the whole thing in time: a piston compresses a fuel–air charge, a spark releases heat near top dead centre, and the hot gas expands — while we watch the pressure, temperature, and flow evolve on a moving mesh.

The simulation solves the compressible Euler equations on a mesh whose top wall is the piston. As the piston moves the cells contract and expand, each cell keeps its total mass/momentum/energy, and the piston boundary applies its p dV work. A Wiebe heat-release law injects the combustion energy near θ = 2π (TDC). From the solved density ρ and pressure p we recover the temperature T = p/(ρ R) (ideal gas), and the flow is the velocity field.

Build the engine and run the cycle

First the geometry: bore B, stroke S, connecting-rod length l, and clearance volume V_c. engine_simulation_2d drives the crank from π (BDC, before compression) through TDC and back to 3π (BDC, after expansion), and records a snapshot of the whole field at regular crank angles:

B, S, l, Vc = 0.08, 0.08, 0.12, 5e-5   # bore, stroke, rod, clearance (m)
γ, R = 1.4, 287.0

θs, Lys, snaps = engine_simulation_2d(30, 30, B, S, l, Vc, 500.0, π, 3π;
                                      γ=γ, R=R, p_bdc=1e5, T_bdc=300.0,
                                      Q=2.0e4, θ_ign=2π, Δθ=0.3, nframes=30)

println(length(snaps), " frames,  θ from ", round(θs[1], digits=2), " to ", round(θs[end], digits=2))

30 frames,  θ from 3.35 to 9.4

θs[k] is the crank angle, Lys[k] the piston height (the piston displacement), and snaps[k] holds the cell-centred fields ρ, u, v, p.

Assemble the animation frames

Each snapshot is a regular 30 × 30 grid of cells. To hand it to the viewer we triangulate each cell into two triangles (a TriangleMesh) and spread each cell’s field value onto its two triangles:

nx, ny = size(snaps[1].p)
Lx = B

function frame_mesh(s, nx, ny, Lx)
    verts = zeros(3, (nx + 1) * (ny + 1))
    for j in 1:ny+1, i in 1:nx+1
        verts[:, i + (nx + 1) * (j - 1)] .= (Lx * (i - 1) / nx, s.Ly * (j - 1) / ny, 0.0)
    end
    faces = zeros(Int, 3, 2 * nx * ny)
    f = 0
    for j in 1:ny, i in 1:nx
        n00 = i + (nx + 1) * (j - 1)
        faces[:, f+1] .= (n00, n00 + 1, n00 + nx + 2)
        faces[:, f+2] .= (n00, n00 + nx + 2, n00 + nx + 1)
        f += 2
    end
    TriangleMesh(verts, faces)
end

function cell_scalar(v, nx, ny)
    out = zeros(2 * nx * ny)
    f = 0
    for j in 1:ny, i in 1:nx
        out[f+1] = v[i, j]; out[f+2] = v[i, j]; f += 2
    end
    out
end

meshes = [frame_mesh(s, nx, ny, Lx) for s in snaps]

30-element Vector{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}}:
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.08937393550531728 0.08937393550531728; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.07080717093762125 0.07080717093762125; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.041575331487259864 0.041575331487259864; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.027593138789830593 0.027593138789830593; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.019975751396124824 0.019975751396124824; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.015006875982048588 0.015006875982048588; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.01265529968998009 0.01265529968998009; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.011223008245369667 0.011223008245369667; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.010286079205909904 0.010286079205909904; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.009986218809332965 0.009986218809332965; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 ⋮
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.025002891083574265 0.025002891083574265; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.030119418825797228 0.030119418825797228; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.03581919491559361 0.03581919491559361; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.043222201694241985 0.043222201694241985; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.05470607955942735 0.05470607955942735; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.06606983764512768 0.06606983764512768; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.07746441498031893 0.07746441498031893; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.0872044353098968 0.0872044353098968; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)
 TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.0026666666666666666 … 0.07733333333333332 0.08; 0.0 0.0 … 0.08997357721883373 0.08997357721883373; 0.0 0.0 … 0.0 0.0], [1 1 … 929 929; 2 33 … 930 961; 33 32 … 961 960], nothing)

Pressure: compression and the combustion spike

vtk(meshes; times=θs, …) plays the frame sequence as an animation. Colour each frame by the pressure to see the charge compress isentropically and then spike when the spark fires:

pfields = [cell_scalar(s.p, nx, ny) for s in snaps]
vtk(meshes; times=θs, field=pfields, fieldname="pressure", colormap=:inferno)

Speed

Zoom Opacity

Temperature

The temperature T = p/(ρ R) rises during compression (isentropic heating) and jumps to thousands of Kelvin as the combustion energy is released:

Tfields = [cell_scalar(s.p ./ (s.ρ .* R), nx, ny) for s in snaps]
vtk(meshes; times=θs, field=Tfields, fieldname="temperature", colormap=:turbo)

Speed

Zoom Opacity

Flow

The gas follows the piston: nearly still at top dead centre, fastest mid-stroke as the piston sweeps through the charge. Colour by the flow speed √(u² + v²):

vfields = [cell_scalar(sqrt.(s.u .^ 2 .+ s.v .^ 2), nx, ny) for s in snaps]
vtk(meshes; times=θs, field=vfields, fieldname="flow speed", colormap=:plasma)

Speed

Zoom Opacity

The piston displacement

The shrinking box height is the piston motion. Lys[k] maps the crank angle to the piston position — down from BDC to TDC, then back out:

for k in (1, 8, 15, 22, length(Lys))
    s = snaps[k]
    println("θ = ", round(θs[k], digits=2),
            "   piston height = ", round(Lys[k] * 1e3, digits=1), " mm",
            "   mean p = ", round(mean(s.p) * 1e-5, digits=1), " bar",
            "   mean T = ", round(mean(s.p ./ (s.ρ .* R)), digits=0), " K")
end

θ = 3.35   piston height = 89.4 mm   mean p = 1.0 bar   mean T = 301.0 K
θ = 6.07   piston height = 11.2 mm   mean p = 18.4 bar   mean T = 690.0 K
θ = 6.58   piston height = 12.3 mm   mean p = 93.0 bar   mean T = 3817.0 K
θ = 7.07   piston height = 25.0 mm   mean p = 34.7 bar   mean T = 2894.0 K
θ = 9.4   piston height = 90.0 mm   mean p = 5.8 bar   mean T = 1734.0 K

The compression stroke (π → 2π) raises the pressure and temperature isentropically; the Wiebe heat release near TDC spikes them to their peak; the expansion stroke (2π → 3π) drops them again as the piston is pushed back down.

Next steps

Internal combustion engine cycle for the 0D cycle, the Otto efficiency, and the moving-mesh compression behind this simulation.
Combustion chemistry for the kinetics and flamelets that close the heat release.
Turbulence closures for the RANS/LES models that feed the turbulent flame speed.