6. Cylinder with Combustion & Valves

The cylinder inherits from the 0‑D volume, adding piston kinematics, valve orifices, and heat release.

Piston Motion

The crank‑slider geometry gives the volume as a function of crank angle \(\theta\):

$$ V(\theta) = V_c + \frac{\pi B^2}{4} \cdot \left[ \frac{S}{2} + l - l\cos\theta - \sqrt{l^2 - \left(\frac{S}{2}\sin\theta\right)^2} \right], $$

where \(B\) is bore, \(S\) stroke, \(l\) rod length, \(V_c\) clearance volume. The first derivative \(dV/dt\) is obtained analytically or via a simple finite difference.

Valve Flow Model

Flow through a poppet valve is modelled as a compressible orifice [1]. The mass flow rate is:

$$ \dot{m} = C_d A_v \frac{p_u}{\sqrt{R_g T_u}} \, \Psi\!\left(\frac{p_d}{p_u}\right), $$

with discharge coefficient \(C_d\), effective area \(A_v\) (function of valve lift), upstream stagnation pressure \(p_u\) and temperature \(T_u\), and the flow function \(\Psi\):

$$ \Psi(r) = \begin{cases} \sqrt{\displaystyle\frac{2\gamma}{\gamma-1}\left(r^{2/\gamma} - r^{(\gamma+1)/\gamma}\right)}, & r > \left(\frac{2}{\gamma+1}\right)^{\gamma/(\gamma-1)},\\ \sqrt{\displaystyle\gamma\left(\frac{2}{\gamma+1}\right)^{(\gamma+1)/(\gamma-1)}}, & \text{choked}. \end{cases} $$

Reverse flow is handled symmetrically. This provides a smooth transition between subsonic and choked conditions.

C# Valve Flow Function

// Pre‑compute constant exponents for performance
double R_gas = 287.058;
double crit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
double choked_psi = Math.Sqrt(gamma * Math.Pow(2.0/(gamma+1.0), (gamma+1.0)/(gamma-1.0)));

double ValveMassFlow(double pUp, double TUp, double pDown, double Cd, double A)
{
    double r = Math.Max(pDown / pUp, 1e-10);
    double psi;
    if (r > crit)
        psi = Math.Sqrt(2.0 * gamma / (gamma - 1.0) *
              (Math.Pow(r, 2.0/gamma) - Math.Pow(r, (gamma+1.0)/gamma)));
    else
        psi = choked_psi;
    return Cd * A * pUp / Math.Sqrt(R_gas * TUp) * psi;
}

Note: In a production build, the Math.Pow calls inside the valve function should be replaced with pre‑computed lookup tables or polynomial approximations for maximum speed.

Combustion

Heat release is modelled using a Wiebe function [2]:

$$ \dot{Q} = Q_{\text{fuel}} \,\frac{a\,(m+1)}{\Delta\theta} \left(\frac{\theta - \theta_0}{\Delta\theta}\right)^m \exp\!\left[-a\left(\frac{\theta - \theta_0}{\Delta\theta}\right)^{m+1}\right], $$

where \(\theta_0\) is spark timing, \(\Delta\theta\) burn duration, \(a\) and \(m\) are shape parameters. For Diesel, an Arrhenius‑based ignition delay triggers the heat release.

C# Wiebe Heat Release

double WiebeHeatRelease(double theta, double theta0, double deltaTheta,
                        double Qfuel, double a, double m)
{
    if (theta < theta0 || theta > theta0 + deltaTheta) return 0.0;
    double x = (theta - theta0) / deltaTheta;
    double rate = (a * (m + 1.0) / deltaTheta) * Math.Pow(x, m) * Math.Exp(-a * Math.Pow(x, m + 1.0));
    return Qfuel * rate;
}

Knock Detection

Knock is detected by monitoring the high‑frequency content of the cylinder pressure signal. The current implementation uses a simple pressure‑derivative threshold as a placeholder; a more complete approach would filter the pressure trace in the knock frequency band (typically 5–15 kHz for automotive engines) and compare the band energy against a calibrated threshold. When knock is flagged, a high‑frequency “ping” sound is triggered via a filtered noise burst, combined with the normal combustion noise.

C# Knock Detection (Placeholder)

bool DetectKnock(double p, double pPrev, double dt, double threshold, double crankAngle, double maxAngle)
{
    double dpdt = (p - pPrev) / dt;
    return (dpdt > threshold) && (crankAngle < maxAngle);
}
// A physically richer indicator would use band‑pass energy of p in the 5–15 kHz range.

References

  1. J. B. Heywood, Internal Combustion Engine Fundamentals, McGraw‑Hill, 2018.
  2. G. P. Blair, Design and Simulation of Four‑Stroke Engines, SAE, 1999.
← PreviousNext →