Skip to content

Recipe: wall temperature, burn shape, and fuelling limits

Three models that each replace a constant you had to guess. Every snippet below runs as written — copy it into a file and execute it.

The theory, and why each is shaped the way it is, is in Theory: wall temperature, burn shape, and fuelling limits.

All three are opt-in

Defaults are unchanged. Turning one on changes results — that is the point — so nothing here is switched on for you.

1. Let the wall find its own temperature

Before: a constant you guessed

from aeolus1d.components.cylinder import Cylinder

cyl = Cylinder(bore=0.130, stroke=0.160, T_wall=450.0)   # a guess

450.0 is right at one operating point and wrong at every other.

After: solve it

from aeolus1d.closures.cylinder_wall_thermal import (
    CylinderWallThermal,
    liner_friction_power,
)

wall = CylinderWallThermal(
    bore=0.130,
    stroke=0.160,
    T_coolant=363.0,      # 90 °C
    T_oil=373.0,          # 100 °C
    tau_cycles=12.0,      # metal lag, in cycles
)

Feed it every solver step, then close the cycle:

# --- once per solver step, inside your march ---
wall.accumulate(
    dt=dt,
    T_gas=T_gas,
    h_head=h, h_piston=h, h_liner=h,     # from your h correlation
)

# --- once per cycle ---
wall.add_friction_power(
    liner_friction_power(
        fmep_pa=1.2e5,
        displacement_m3=2.12e-3,
        cycles_per_s=rpm / 120.0,
    )
)
surfaces = wall.update()
print(surfaces.T_head, surfaces.T_piston, surfaces.T_liner)

update() returns the surfaces for the next cycle and resets the integrals. Hand them to the cylinder's zone temperatures.

A complete, runnable demonstration

import math
from aeolus1d.closures.cylinder_wall_thermal import CylinderWallThermal

def run(T_gas_peak: float, n_cycles: int = 400) -> float:
    wall = CylinderWallThermal(bore=0.130, stroke=0.160, tau_cycles=12.0)
    for _ in range(n_cycles):
        for i in range(720):                      # 1° steps
            theta = math.radians(i)
            # crude stand-ins: hot and high-h near TDC
            hot = math.exp(-((i - 370) / 40.0) ** 2)
            wall.accumulate(
                dt=1.0e-5,
                T_gas=400.0 + T_gas_peak * hot,
                h_head=200.0 + 2500.0 * hot,
                h_piston=200.0 + 2500.0 * hot,
                h_liner=150.0 + 1500.0 * hot,
            )
        wall.update()
    return wall.surfaces.T_head

for peak in (1200.0, 1800.0, 2400.0):
    print(f"peak gas {peak:6.0f} K  ->  head {run(peak):6.1f} K")

The head temperature rises with load, which a constant cannot do. That is the whole point of the model.

Start with friction alone

Before wiring the full network, test the cheap hypothesis:

from aeolus1d.closures.cylinder_wall_thermal import liner_friction_power

q = liner_friction_power(fmep_pa=1.2e5, displacement_m3=2.12e-3,
                         cycles_per_s=12.5)
print(f"{q / 1000.0:.2f} kW into the liner, per cylinder")

If adding that to your existing constant-wall closure does not move a Q_WHT shortfall in the right direction and with the right load trend, the resistance network will not rescue it either. Cheaper to find out now.

fraction_to_liner=0.25 is a starting value

It is the older literature's ring-friction split, unverified here. Measure it before trusting it.

2. Make the burn shape follow the ignition delay

One field

from aeolus1d.closures.diesel_predictive_combustion import (
    DieselPredictiveCombustion,
)

model = DieselPredictiveCombustion(
    soi_deg=346.0,
    inj_duration_deg=18.0,
    m_fuel=1.2e-4,
    burn_law="watson",          # <- was "wiebe"
)

That is the whole change. burn_law="wiebe" (the default) keeps a fixed premixed duration; "watson" evaluates SAE 800029 directly, so the premixed sharpness follows τ_id · N and the diffusion rate follows φ.

Seeing it respond

import math
from aeolus1d.closures.diesel_predictive_combustion import (
    DieselPredictiveCombustion,
)

omega = 2.0 * math.pi * 1400.0 / 60.0
common = dict(m_charge=2.7e-3, omega_crank=omega, mean_piston_speed=8.0)

model = DieselPredictiveCombustion(soi_deg=346.0, burn_law="watson")

for cycle, (p, T, label) in enumerate(
    [(90e5, 1000.0, "hot charge"), (35e5, 680.0, "cold charge")], start=1
):
    model.advance(1e-5, 360.0, cycle=cycle, p=p, T=T, **common)
    print(f"{label:12s} delay {model.last_tau_id_deg:5.2f}°CA  "
          f"K1 {model._msw.shape.K1:.3f}  beta {model.last_premixed_fraction:.3f}")

The colder charge delays ignition, and K1 follows it. A fixed premixed duration cannot express that, and phases wrongly as soon as load moves the delay.

Using the law on its own

from aeolus1d.closures.watson_burn_law import (
    watson_burn_rate, watson_burnt_fraction, watson_shape,
)

shape = watson_shape(tau_id_ms=1.0, rpm=1500.0, phi=0.5)
for tau in (0.0, 0.02, 0.05, 0.2, 0.5, 1.0):
    print(f"tau {tau:4.2f}  x_b {watson_burnt_fraction(tau, 0.4, shape):.4f}"
          f"  dx/dtau {watson_burn_rate(tau, 0.4, shape):8.3f}")

tau_id_ms is milliseconds

Passing seconds does not raise — it collapses K1 to 2.0 and silently disables the correlation. The difference is only about 1.5×, which is why it is easy to miss.

When the injection rate matters, use Chmela instead

Watson cannot see the injector. A pilot, a post and a ramped needle all give it the same curve, because the correlation's only inputs are the delay and φ. burn_law="chmela" reads an actual injection profile:

from aeolus1d.closures.chmela_mcc import InjectionRate
from aeolus1d.closures.diesel_predictive_combustion import (
    DieselPredictiveCombustion,
)

dps = 1400.0 * 6.0          # deg/s at 1400 rpm
m_fuel = 1.2e-4

pilot_and_main = InjectionRate(
    theta_deg=(339.9, 340.0, 343.0, 343.1, 345.9, 346.0, 364.0, 364.1),
    mdot=(0.0, 0.25 * m_fuel * dps / 3.0, 0.25 * m_fuel * dps / 3.0, 0.0,
          0.0, 0.75 * m_fuel * dps / 18.0, 0.75 * m_fuel * dps / 18.0, 0.0),
)

model = DieselPredictiveCombustion(
    soi_deg=346.0, inj_duration_deg=18.0, m_fuel=m_fuel,
    burn_law="chmela",
    injection_rate=pilot_and_main,
)

Without injection_rate it falls back to a rectangle built from soi_deg/inj_duration_deg — it runs, but a rectangle has no rate shape, so it gives up most of what Chmela offers over Watson.

Seeing it respond:

import math
from aeolus1d.closures.chmela_mcc import ChmelaMCC, rectangular_injection

OMEGA = 2.0 * math.pi * 1500.0 / 60.0
M = 1.2e-4

def burn(duration_deg: float, hole_m: float = 1.6e-4) -> tuple[float, float]:
    inj = rectangular_injection(soi_deg=350.0, duration_deg=duration_deg,
                                m_fuel=M, omega_crank=OMEGA)
    m = ChmelaMCC(injection=inj, m_fuel=M, hole_diameter_m=hole_m)
    m.reset_cycle()
    dps, n = OMEGA * 180.0 / math.pi, 4000
    dth = 140.0 / n
    dt, theta, peak = dth / dps, 340.0, 0.0
    for _ in range(n):
        V = 8e-5 + 9e-4 * (1.0 - math.cos(math.radians(theta - 360.0)))
        peak = max(peak, m.advance(dt, theta, V=max(V, 8e-5),
                                   omega_crank=OMEGA))
        theta += dth
    return m.burnt_fraction, peak / 1e3

for label, kw in (("20 deg, 0.16 mm", {}),
                  ("10 deg, 0.16 mm", {"duration_deg": 10.0}),
                  ("20 deg, 0.12 mm", {"hole_m": 1.2e-4})):
    f, pk = burn(kw.pop("duration_deg", 20.0), **kw)
    print(f"{label:18s} burnt {f:.3f}  peak {pk:7.1f} kW")

A shorter injection and a smaller nozzle both raise the peak, through the spray turbulence they generate. Nothing in Watson can express either.

The constants are calibration, and two were not verified

C_mod and C_rate are calibration constants in every published use of this model. C_turb and C_diss govern the turbulence sub-model and were not verified against the source paper — they are exposed as parameters rather than buried for exactly that reason. Calibrate against measured heat release.

It integrates; the others sample

Chmela carries released heat and spray kinetic energy across the cycle. Sampling the same crank angle twice does not give the same answer, as it would for a Wiebe or Watson curve. reset_cycle() starts a fresh cycle; the closure does that for you on a new cycle index.

3. Cap fuelling at what the air can burn

The chain

import math
from aeolus1d.control import SpeedGovernor, SmokeLimiter, governed_fuelling

governor = SpeedGovernor(
    m_fuel_max=2.0e-4,                          # kg/cycle/cylinder at the stop
    omega_set=2.0 * math.pi * 1500.0 / 60.0,
    droop=0.05,                                 # 5 % droop
)
limiter = SmokeLimiter(lambda_min=1.4, afr_stoich=14.5)

Then once per cycle:

result = governed_fuelling(
    governor, limiter,
    omega=omega_now,
    dt=cycle_time_s,
    m_air_trapped=m_charge * (1.0 - burnt_fraction),   # AIR, not charge
)
print(f"demand {result.m_fuel_demand:.3e}  ->  fuelled {result.m_fuel:.3e}"
      f"  limited={result.limited}  lambda={result.lambda_actual:.2f}")

A load step, where it earns its place

import math
from aeolus1d.control import SpeedGovernor, SmokeLimiter, governed_fuelling

omega_set = 2.0 * math.pi * 1500.0 / 60.0
governor = SpeedGovernor(m_fuel_max=2.0e-4, omega_set=omega_set, droop=0.05)
limiter = SmokeLimiter(lambda_min=1.4)

# the turbo spools over ~1 s; trapped air follows it
for i in range(6):
    m_air = 1.0e-3 + i * 1.4e-3
    r = governed_fuelling(governor, limiter, omega=omega_set * 0.93,
                          dt=0.04, m_air_trapped=m_air)
    flag = "LIMITED" if r.limited else "free"
    print(f"air {m_air * 1e3:5.2f} g  fuel {r.m_fuel * 1e6:6.1f} mg  "
          f"lambda {r.lambda_actual:5.2f}  {flag}")

Early on the governor wants full fuel and the air refuses; as boost arrives the cap lifts. Without the limiter every one of those cycles would burn fuel that has no oxygen.

Feed it air, not charge

Recirculated exhaust is mass and is not oxygen. Passing total charge mass at a high-EGR point licenses exactly the over-fuelling the limiter exists to prevent — worst where it matters most.

4. Share a component between cases

A part lives in its own file with no [case] block. Bundled parts live under examples/parts/, kept out of the case lineup on purpose: a part has no [case] block and cannot be run, so aeolus1d examples does not list one and nothing counts it as a bundled case.

# examples/parts/std_pipe.toml
length  = 0.75
n_cells = 200
area    = 0.0035

[init]
p = 1.0e5
T = 300.0
u = 0.0

A case names it and overrides only what differs:

[[pipe]]
id     = "tube"
ref    = "parts/std_pipe"
length = 1.25          # n_cells, area and [init] come from the part

from aeolus1d.io.parsers import load_case

case = load_case("examples/92_component_library.toml")
print(case.pipes[0].length)      # 1.25  — the override
print(case.pipes[0].n_cells)     # 200   — from the part
print(case.component_refs)       # {'pipe[0]': 'parts/std_pipe'}

Local keys always win. The part supplies defaults; the case overrides them. The reverse would let a re-measured part silently change a case that had deliberately overridden it.

Refs resolve against the case directory first, then library_roots, then $AEOLUS_COMPONENT_LIBRARY (a PATH-style list) — so a case can ship a local override of a shared part without touching the library:

export AEOLUS_COMPONENT_LIBRARY=/srv/parts:/home/me/parts

Nested tables merge; arrays are replaced wholesale. A part may itself carry ref, so a part can be a variant of another part. Cycles are refused with the chain named.

Saving flattens the reference

The GUI writer emits every field explicitly, so saving a case built from a ref inlines the part and loses the link. It still reloads to an identical case, but the next re-measurement of that part will no longer reach it. save_case warns when this happens — keep the original TOML as the source of truth.

Run the bundled example:

aeolus1d run examples/92_component_library.toml

Where to go next