Theory: wall temperature, burn shape, and fuelling limits¶
Three models that share a theme — each replaces a constant the user had to guess with something the simulation computes from state it already has.
| Model | Was | Is |
|---|---|---|
| Wall temperature | one constant T_wall |
head / piston / liner, solved per cycle |
| Burn shape | fixed premixed duration | shape follows the ignition delay |
| Fuelling | prescribed m_fuel |
speed → demand → λ cap |
They are all opt-in. Every number recorded on this project was produced without them, and switching a default silently would invalidate all of it at once.
1. Cylinder wall temperature¶
The problem¶
Cylinder takes T_wall — optionally split into T_head, T_piston,
T_liner — as an input. A real surface temperature is an output of the
heat the gas puts into it, and moves 50–150 K across a load range.
The symptom on the EPW diesel: the model under-predicts wall heat Q_WHT at
31.8 kW against a measured 46.7 kW at 1200 rpm, and — the part that matters —
the shortfall collapses with load:
| op | shortfall, as % of fuel power |
|---|---|
| 1200 rpm | 6.45 % |
| 1500 rpm | 1.87 % |
| 1800 rpm | 0.78 % |
A wall that cannot warm up cannot produce a load-dependent shortfall in either direction. Tuning the constant fixes one operating point and breaks the others.
The model¶
Three surface nodes, each behind a gas-side convective resistance, joined through metal to a coolant or oil sink:
T_gas
│ R_gas = 1/(h̄·A) rebuilt every cycle
┌────┴────┐
│ surface │ ← friction, on the liner only
└────┬────┘
│ R_metal = t/(k·A)
│ R_sink = 1/(h_sink·A)
T_coolant (head, liner) / T_oil (piston)
Solving the node balance for each surface:
with \(R_\text{back} = R_\text{metal} + R_\text{sink}\).
This is the formulation of Torregrosa, Olmeda, Degraeuwe & Reyes (2006).
Its contribution over the older textbook network is that every resistance is
written as a function of geometry, operating parameters and material
properties, so it transfers between engines rather than being re-fitted for
each one. Nothing in CylinderWallThermal is a bare number; each resistance
is built from bore, stroke, thickness and conductivity.
Three details that are easy to get wrong¶
The gas temperature is h-weighted, not time-weighted. Both \(\int h\,dt\) and \(\int h\,T_\text{gas}\,dt\) are accumulated, giving
The gas is hottest exactly when \(h\) is largest. A plain time average under-weights the part of the cycle doing nearly all the heating, and reports a cooler wall than the cycle actually produces.
The liner area is half the stroke. The liner is fully exposed at BDC and barely at TDC, while \(h\) peaks near TDC. Half-stroke is the cycle-average that pairs with a cycle-averaged coefficient; the full stroke roughly doubles the liner's share of the heat.
The metal lags. Its time constant is seconds; a cycle is milliseconds. The
model solves the steady network each cycle and relaxes toward it over
tau_cycles. Jumping straight to the steady answer would make the readout
depend on which cycle it was taken at.
Friction, kept separable on purpose¶
A share of piston-ring friction power heats the liner rather than leaving with the oil:
liner_friction_power is a free function, not a method, so it can be
tried against the existing constant-wall closure on its own. That is the
cheapest possible test of whether the missing Q_WHT is friction at all: if
adding this term does not move the shortfall in the right direction and with
the right load trend, the resistance network will not rescue it either.
0.25 is a starting value, not a measurement
DEFAULT_FRICTION_TO_LINER = 0.25 is the older literature's ring-friction
split. It has not been verified on this project. Measure it before
relying on it.
Choosing a gas-side correlation¶
| Model | Length scale | T exponent | Notes |
|---|---|---|---|
WoschniHeatTransfer |
bore (constant) | −0.55 | the default; C2 term opt-in |
HohenbergHeatTransfer |
volume | −0.40 | cheaper, degenerates gracefully at zero piston speed |
ChangHeatTransfer |
instantaneous chamber height | −0.73 | combustion term ÷ 6 |
Chang et al. (2004) was derived from instantaneous surface heat-flux measurements rather than a cycle-averaged energy balance, which is why its temperature exponent disagrees with Woschni's. Its length scale \(L(t) = V/(\pi B^2/4)\) thins as the piston approaches TDC, where Woschni's \(B^{-0.2}\) is constant over the cycle. And it divides Woschni's combustion-turbulence term by 6 — unmodified, that term over-predicts the peak against measured flux.
2. The Watson burn law¶
The problem¶
The predictive closure already had Watson's β — how much of the fuel burns
premixed. What it could not predict was the shape: premixed_duration_deg
was fixed at 8°, and the two Wiebe stages had constant exponents.
A premixed spike whose width cannot respond to the ignition delay phases wrongly the moment load moves the delay. That is the mechanism behind the recorded defect "MFB50 runs 2.4–2.7 CAD early at mid and high load."
The correlation¶
Watson, Pilley & Marzouk, SAE 800029 (1980), with \(\tau = (\theta - \theta_\text{soc})/\Delta\theta_\text{burn} \in [0,1]\):
\(K_1\) is the point. It rises with \(\tau_\text{id} N\) — the delay expressed in crank angle — so a longer delay gives a later, sharper premixed spike. \(K_3\) falls with φ, so a richer charge burns its diffusion phase faster.
\((\tau_\text{id} N)\) is in ms·rpm
Passing seconds does not raise. It collapses \(K_1\) to its constant term 2.0, which looks like a plausible exponent while having silently disabled the correlation. The ratio is only about 1.5 — which is exactly why the unit error is easy to miss by eye.
It is not a Wiebe, and cannot be made into one¶
A Wiebe with the usual completeness constant \(a = 6.908\) reaches half its mass at \(\tau_{1/2} = (\ln 2/a)^{1/(m+1)}\), which for any \(m \geq 0\) is at least 0.100. Watson's premixed spike at a representative heavy-duty point (\(K_1 = 2.52\), \(K_2 = 5000\)) reaches half its mass at 0.030. Matching it needs \(m = -0.35\) — an infinite burn rate at ignition.
A "Wiebe-equivalent exponent" helper was written for this project, measured, found ill-posed rather than merely inaccurate, and deleted. Callers wanting Watson's phasing must evaluate Watson's own law.
Numerics¶
\(K_2 = 5000\) makes \((1 - \tau^{K_1})^{K_2}\) underflow to 0.0 — or raise for \(\tau^{K_1} > 1\) — well inside the range the model is evaluated over. Both the fraction and the rate go through \(\exp(K_2 \ln(1-x))\) instead.
The burn rate is analytic, not a finite difference:
Differencing two underflowed fractions returns noise or exactly zero, which erases the spike — the entire reason the correlation is here.
2b. Chmela: mixing-controlled heat release¶
Watson is a correlation: given a delay and a φ it produces a burn curve. Nothing in it knows the injector did anything after SOI, so a pilot, a post, a rate shape and a rail-pressure change all produce the same curve. On a common-rail engine that is most of what the calibration actually moves.
Chmela & Orthaber, SAE 1999-01-0186 (1999) asks a different question. For diffusion combustion the rate is not set by chemistry — it is set by how fast fuel and air mix, and two things control that, both from the injection:
\(f_1\) is in kg and \(f_2\) in 1/s, so the product is a fuel burn rate and the explicit \(H_u\) makes it watts. Sources fold \(H_u\) in differently; it is written out so the units check by eye.
The turbulence comes from the spray momentum itself:
What is verified, and what is not¶
Read this before calibrating
The structure above is well established and reproduced consistently across the literature — the two controlling quantities, \(C_\text{mod}\) and \(C_\text{rate}\), and the \(\sqrt{k}/V^{1/3}\) mixing frequency.
The numerical constants are a different matter. \(C_\text{mod}\) and \(C_\text{rate}\) are calibration constants in every published use of this model — that is how it is meant to be used. \(C_\text{turb}\) and \(C_\text{diss}\) govern the turbulence sub-model, and the values defaulted here were not verified against the source paper. They are exposed and documented as such rather than buried, because a fitted constant that reads like a published one is how a model acquires false authority.
C_mod = 0.30 was chosen by sweeping a nominal heavy-duty point — 1500 rpm,
1.2e-4 kg, a 20° rectangular injection from 350°, an 8 × 0.16 mm nozzle, a
340–480° window:
C_mod |
burnt | C_mod |
burnt |
|---|---|---|---|
| 0.03 | 0.264 | 0.12 | 0.705 |
| 0.05 | 0.399 | 0.20 | 0.868 |
| 0.10 | 0.638 | 0.30 | 0.952 |
A smaller default is not conservative, it is broken: at 1e-3 the model releases a fraction of a percent and reads as a bug rather than a calibration task.
Diffusion only¶
Chmela's model deliberately does not describe the premixed spike — that phase
is kinetically controlled, not mixing-controlled. ChmelaMCC releases the
premixed share through a separate fast Wiebe so the total is conserved.
Releasing it through the mixing law would make it respond to spray
turbulence, which is the wrong physics for the one phase that is chemistry.
It is a rate model, not a curve¶
Watson and the Wiebe evaluate \(x_b(\theta)\) — sampling the same crank angle twice gives the same answer. Chmela integrates, carrying released heat and spray kinetic energy across the cycle.
That difference produced a real bug, found by a test and worth recording. When a new cycle began but SOI had not yet been reached, the previous cycle's burn was still being advanced. For a curve that is invisible — evaluated before its own SOC it returns ~0 either way. For a rate model it is not: the spent burn still held \(f_1 = 1.0\times10^{-6}\) kg of unburnt fuel and a live spray kinetic energy, and dribbled 42.6 J into the following cycle — 0.8 % of its energy, at the wrong crank angle. The closure now drops the stale burn on a new cycle index.
Choosing between the three¶
| responds to delay | responds to φ | responds to injection rate | |
|---|---|---|---|
"wiebe" |
SOC only | β only | ✗ |
"watson" |
✓ shape | ✓ shape | ✗ |
"chmela" |
via β | via β | ✓ |
Watson is cheap and needs only delay and φ — use it when that is all the case carries, and it is what the recorded MFB50 defect calls for. Chmela needs an injection-rate profile; without one it falls back to a rectangle, which runs but gives up most of the benefit.
3. Governor and smoke limiter¶
The problem¶
Fuelling is prescribed: m_fuel per cycle is an input. That is correct for a
steady operating-point sweep and wrong for anything transient, because it
leaves out the loop that decides how much fuel an engine actually burns.
Not a rack governor¶
The 1980s form of this is flyweight dynamics driving a normalised rack-position pump map. Production diesels have not worked that way in decades, and a rack position is an actuator coordinate with no meaning in a model whose injector is commanded in kg. What replaced it:
ω → SpeedGovernor → demand → SmokeLimiter → m_fuel
↑
m_air_trapped
SpeedGovernor is droop or isochronous. droop is a machine
specification, not a gain: 5 % droop means full fuel at 5 % under target.
droop = 0 is isochronous — pure integral, holds exact speed — and the
integrator is clamped, because a missing anti-windup shows up only after the
error reverses, as an engine that will not come off full fuel.
SmokeLimiter is the half that matters:
This is the standard ECU function. Without it a transient fuel demand outruns the air system — the turbocharger cannot spool instantly — and the model burns fuel that has no oxygen, producing torque the real engine would refuse and smoke it would emit.
Three decisions with silent failure modes¶
λ, not AFR. λ is dimensionless and transfers across fuels; AFR does not. A limiter written in AFR mis-calibrates the moment the fuel changes.
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.
Callers holding a charge mass and a burnt fraction should pass
m_charge · (1 − burnt_fraction).
lambda_min < 1.0 raises. Below stoichiometric there is not enough oxygen
to burn the fuel at all, and a "limit" there is not a limit.
Which constraint bound?¶
FuellingResult reports limited. Without it a torque shortfall caused by
the limiter and one caused by the governor are indistinguishable in the
output, and the wrong one gets tuned.
The order is not interchangeable: the governor asks what the speed error wants, the limiter refuses what the air cannot support. Reversed, the integrator winds up against a cap it cannot see and the engine over-fuels the moment air arrives.
References¶
- A. Torregrosa, P. Olmeda, B. Degraeuwe, M. Reyes, "A concise wall temperature model for DI Diesel engines", Applied Thermal Engineering 26 (2006) 1320–1327.
- Chang et al., "New Heat Transfer Correlation for an HCCI Engine Derived from Measurements of Instantaneous Surface Heat Flux", SAE 2004-01-2996.
- G. Woschni, SAE 670931 (1967). · G. Hohenberg, SAE 790825 (1979).
- N. Watson, A. D. Pilley, M. Marzouk, "A Combustion Correlation for Diesel Engine Simulation", SAE 800029 (1980).
- N. Watson, M. S. Janota, Turbocharging the Internal Combustion Engine, Macmillan 1982, §15.
- F. Chmela, G. Orthaber, SAE 1999-01-0186. · M. Barba et al., SAE 2000-01-2933.