Internal Combustion Engine Cycle: Kinematics, Thermodynamics, and Combustion

Overview

An internal combustion engine converts the heat of burning fuel into work. In the four-stroke spark-ignition cycle a piston slides up and down a cylinder: it compresses a fuel–air charge, a spark ignites it, the hot gases expand and push the piston, and the burnt charge is exhausted. Rasmah models the 0D backbone of this process — the slider-crank piston motion, the cylinder volume, the ideal Otto cycle, and the Wiebe burn law — then assembles them into a firing in-cylinder cycle whose efficiency is checked against the classical formula.

Slider-crank kinematics

The piston is driven by a crank of radius a (half the stroke) and a connecting rod of length l. Its position is x(θ) = a cos θ + √(l² − a² sin² θ), so at TDC (θ = 0) the piston is at the top and at BDC (θ = π) at the bottom:

a, l = 0.04, 0.12                    # crank radius, rod length (m)
println("displacement at TDC = ", piston_displacement(0.0, a, l), " m  (0)")
println("displacement at BDC = ", piston_displacement(π, a, l), " m  (= stroke ", 2a, ")")

displacement at TDC = 0.0 m  (0)
displacement at BDC = 0.08000000000000002 m  (= stroke 0.08)

The instantaneous cylinder volume adds the swept volume of the piston to the clearance volume V_c:

B, S, Vc = 0.08, 0.08, 5e-5          # bore, stroke, clearance volume
V_tdc = cylinder_volume(0.0, B, S, l, Vc)
V_bdc = cylinder_volume(π, B, S, l, Vc)
r = compression_ratio(B, S, Vc)
println("V_TDC = ", V_tdc, "  V_BDC = ", round(V_bdc, digits=8), "  r = ", round(r, digits=3))

V_TDC = 5.0e-5  V_BDC = 0.00045212  r = 9.042

The ideal Otto cycle

The ideal spark-ignition cycle compresses isentropically, adds heat at constant volume (the spark), expands isentropically, and rejects heat. Its thermal efficiency is η = 1 − 1/r^{γ−1}:

for r in (8.0, 10.0, 12.0)
    println("r = ", r, "   η = ", round(otto_efficiency(r, 1.4), digits=4))
end
# higher compression ratio → higher efficiency (the textbook result)

r = 8.0   η = 0.5647
r = 10.0   η = 0.6019
r = 12.0   η = 0.6299

The Wiebe burn law

Real combustion is not instantaneous — the fuel burns over a finite crank angle. The Wiebe law gives the mass fraction burned:

x_b(θ) = 1 − exp(−a ((θ − θ₀)/Δθ)^{m+1}),

rising from 0 at ignition θ₀ to essentially 1 over the burn duration Δθ:

println("x_b at θ₀      = ", wiebe_burn_fraction(0.0, 0.0, 10.0))
println("x_b at θ₀ + Δθ = ", round(wiebe_burn_fraction(10.0, 0.0, 10.0), digits=4), "  (= 1 − e⁻⁵)")

x_b at θ₀      = 0.0
x_b at θ₀ + Δθ = 0.9933  (= 1 − e⁻⁵)

The firing in-cylinder cycle

Assembling these, single_zone_cycle integrates the single-zone energy equation m c_v dT/dθ = −p dV/dθ + Q dx_b/dθ through compression, Wiebe combustion, and expansion, giving the pressure trace. The area of the P–V diagram is the indicated work, and η = W/Q is the thermal efficiency:

θs = range(π, 3π, length=4000)                       # BDC → TDC → BDC
γ, R = 1.4, 287.0
p_bdc, T_bdc = 1e5, 300.0

_, p_motored = single_zone_cycle(θs, B, S, l, Vc, p_bdc, T_bdc, γ, R, 0.0, 2π, 1.0)
println("motored peak p = ", round(maximum(p_motored) * 1e-6, digits=2), " MPa  (= p_BDC r^γ = ",
        round(motored_peak_pressure(p_bdc, r, γ) * 1e-6, digits=2), " MPa)")

_, p_fired = single_zone_cycle(θs, B, S, l, Vc, p_bdc, T_bdc, γ, R, 750.0, 2π, 0.5)
W = indicated_work(θs, p_fired, B, S, l, Vc)
η = cycle_efficiency(θs, p_fired, B, S, l, Vc, 750.0)
println("fired peak p  = ", round(maximum(p_fired) * 1e-6, digits=2), " MPa   (combustion raises it)")
println("indicated work = ", round(W, digits=1), " J    η = ", round(η, digits=4), "  (Otto limit ", round(otto_efficiency(r, γ), digits=4), ")")

motored peak p = 2.18 MPa  (= p_BDC r^γ = 2.18 MPa)
fired peak p  = 4.81 MPa   (combustion raises it)
indicated work = 409.7 J    η = 0.5463  (Otto limit 0.5855)

The firing cycle does net positive work and its efficiency approaches the ideal Otto value as the burn duration shortens toward the constant-volume limit.

Crank-driven compression (moving-mesh Euler)

The 0D cycle above tracks only the volume. To resolve the gas in space, Rasmah solves the compressible Euler equations on a mesh that moves with the piston: the mesh contracts uniformly behind the moving wall, the conservative update advances each cell’s total mass/momentum/energy, and the piston boundary applies its p dV work. For a slowly moving piston the compression is isentropic, so the pressure obeys p = p₀ (L₀/L)^γ:

r1d = piston_compression_1d(200, 0.02, 20.0)          # nx, piston speed, time
println("mass  = ", sum(r1d.ρ) * r1d.L / 200, "  (= ρ₀ L₀ = 1)")
println("p     = ", round(sum(r1d.p) / 200, digits=4), "   isentropic p = ", round((1 / r1d.L)^1.4, digits=4))

mass  = 0.9999999999999997  (= ρ₀ L₀ = 1)
p     = 2.0446   isentropic p = 2.0445

Finally, motored_compression_1d drives that same moving-mesh solve with the actual slider-crank piston speed v_p = −ω dL/dθ through the compression stroke, and its peak (TDC) pressure recovers the 0D motored value p_BDC r^γ:

θc, pc = motored_compression_1d(200, B, S, l, Vc, 1000.0, π, 2π; γ=1.4, p_bdc=1e5)
println("peak p = ", round(pc[end] * 1e-6, digits=2), " MPa   isentropic = ",
        round(motored_peak_pressure(1e5, r, 1.4) * 1e-6, digits=2), " MPa")

peak p = 2.18 MPa   isentropic = 2.18 MPa

The firing cycle in space (2D/3D moving mesh)

engine_cycle_2d (and engine_cycle_3d) drive the moving-mesh compressible solver with the slider-crank piston through the whole compression → combustion → expansion stroke, injecting the Wiebe heat release near TDC. The motored run is reversible (the pressure returns to p_bdc and the work vanishes); the fired run does net work with an efficiency approaching the Otto limit:

θm, pm, Wm = engine_cycle_2d(16, 16, B, S, l, Vc, 500.0, π, 3π; γ=1.4, p_bdc=1e5)
println("motored: W = ", round(Wm, digits=2), " J (≈0),  peak = ",
        round(maximum(pm) * 1e-6, digits=2), " MPa")

Q = 2.0e4
θf, pf, Wf = engine_cycle_2d(16, 16, B, S, l, Vc, 500.0, π, 3π; γ=1.4, p_bdc=1e5, Q=Q, θ_ign=2π, Δθ=0.3)
println("fired:   W = ", round(Wf, digits=1), " J,  η = ", round(Wf / Q, digits=3),
        "  (Otto ", round(otto_efficiency(r, 1.4), digits=3), ")")

motored: W = 1.94 J (≈0),  peak = 2.17 MPa
fired:   W = 11394.5 J,  η = 0.57  (Otto 0.586)

The full four-stroke cycle (valves)

A true four-stroke engine also exchanges gas. four_stroke_1d drives the same moving-mesh solver through two crank revolutions, opening a pressure-outlet valve in the head for the intake (p_intake) and exhaust (p_exhaust) strokes and closing it for compression and power. The charge is drawn in, compressed, burned, expanded, and vented — the complete gas-exchange loop:

θ4, p4, W4 = four_stroke_1d(200, B, S, l, Vc, 300.0; γ=1.4, R=287.0, Q=2.0e4, θ_ign=2π, Δθ=0.3)
println("intake p  = ", round(p4[argmin(abs.(θ4 .- π))] * 1e-5, digits=2), " bar  (= p_intake)")
println("peak  p   = ", round(maximum(p4) * 1e-5, digits=1), " bar   (isentropic ",
        round(1e5 * r^1.4 * 1e-5, digits=1), " bar + combustion)")
println("exhaust p = ", round(p4[argmin(abs.(θ4 .- 4π))] * 1e-5, digits=2), " bar  (= p_exhaust)")
println("η = ", round(W4 / 2.0e4, digits=3), "  (Otto ", round(otto_efficiency(r, 1.4), digits=3), ")")

intake p  = 1.02 bar  (= p_intake)
peak  p   = 24.3 bar   (isentropic 21.8 bar + combustion)
exhaust p = 0.97 bar  (= p_exhaust)
η = 0.544  (Otto 0.586)

Next steps

• Combustion chemistry for the Arrhenius kinetics, NASA thermodynamics, and flamelets that close the burn rate.
• Turbulence closures for the RANS/LES models that provide the turbulence intensity feeding the turbulent flame speed.