crank case based two stroke testing

Components/Crankcase.cs

@@ -0,0 +1,150 @@
+// ============================================================
+// File: Crankcase.cs
+// ============================================================
+using System.Collections.Generic;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    public class Crankcase : IComponent
+    {
+        private readonly Crankshaft _crankshaft;
+        private readonly float _crankRadius, _conrodLength, _pistonArea;
+        private readonly float _clearanceVolume, _obliquity;
+
+        private float _mass, _internalEnergy, _airFraction;
+        public float Pressure { get; private set; }
+        public float Temperature { get; private set; }
+        public float Density => _mass / MathF.Max(Volume, 1e-12f);
+        public float Volume { get; private set; }
+
+        public Port IntakePort { get; }
+        public Port TransferPort { get; }
+        private readonly List<Port> _ports;
+        public IReadOnlyList<Port> Ports => _ports;
+
+        // FIX: store previous volume to calculate PdV work
+        private float _prevVolume;
+
+        private const float Rgas = 287.0f;
+        private const float Gamma = 1.4f;
+        private const float Cv = Rgas / (Gamma - 1.0f);
+
+        public Crankcase(Crankshaft crankshaft,
+                         float crankRadius, float conrodLength, float bore,
+                         float clearanceVolume,
+                         float initialPressure, float initialTemperature)
+        {
+            _crankshaft = crankshaft;
+            _crankRadius = crankRadius;
+            _conrodLength = conrodLength;
+            _pistonArea = MathF.PI * 0.25f * bore * bore;
+            _clearanceVolume = clearanceVolume;
+            _obliquity = crankRadius / conrodLength;
+
+            Pressure = initialPressure;
+            Temperature = initialTemperature;
+            float rho = initialPressure / (Rgas * initialTemperature);
+            _mass = rho * clearanceVolume;
+            _internalEnergy = _mass * Cv * initialTemperature;
+            _airFraction = 1.0f;
+            Volume = clearanceVolume;
+            _prevVolume = Volume;
+
+            IntakePort = new Port { Owner = this };
+            TransferPort = new Port { Owner = this };
+            _ports = new List<Port> { IntakePort, TransferPort };
+        }
+
+        public void PreStep(float dt)
+        {
+            // Save previous volume before updating
+            _prevVolume = Volume;
+
+            float theta = _crankshaft.CrankAngleRad % (2f * MathF.PI);
+            float cosTh = MathF.Cos(theta);
+            float sinTh = MathF.Sin(theta);
+            float term = MathF.Sqrt(1f - _obliquity * _obliquity * sinTh * sinTh);
+
+            // FIX: correct piston displacement: downstroke reduces crankcase volume
+            float x = _crankRadius * (1f - cosTh) + _conrodLength * (1f - term);
+            // Maximum volume at TDC (x = 0), minimum at BDC (x = stroke)
+            float maxVolume = _clearanceVolume + _pistonArea * 2f * _crankRadius; // stroke = 2 * crankRadius
+            Volume = maxVolume - _pistonArea * x;
+
+            // Update thermodynamic state using the new volume (before mass transfer)
+            if (_mass > 1e-12f && Volume > 1e-12f)
+            {
+                Temperature = _internalEnergy / (_mass * Cv);
+                Pressure = _mass * Rgas * Temperature / Volume;
+            }
+        }
+
+        public void UpdateState(float dt)
+        {
+            // ---- Mass and energy transport (identical to original) ----
+            float mdotIn = IntakePort.MassFlowRate;
+            float mdotOut = TransferPort.MassFlowRate;
+
+            float dm = (mdotIn - mdotOut) * dt;
+            float dE = (mdotIn * IntakePort.SpecificEnthalpy
+                      - mdotOut * (Cv * Temperature + Pressure / MathF.Max(Density, 1e-12f))) * dt;
+            float dY = (mdotIn * IntakePort.AirFraction - mdotOut * _airFraction) * dt;
+
+            _mass += dm;
+            _internalEnergy += dE;
+            if (_mass > 1e-12f)
+                _airFraction = Math.Clamp((_airFraction * (_mass - dm) + dY) / _mass, 0f, 1f);
+
+            // ---- FIX: add mechanical work done BY the gas ON the piston ----
+            // During a step, volume changed from _prevVolume to current Volume.
+            // Work done BY gas = P * dV (if dV > 0, gas expands and does work, losing energy)
+            float dV = Volume - _prevVolume;
+            // Use average pressure during the step (approximate with current pressure)
+            _internalEnergy -= Pressure * dV;   // removes energy when volume increases
+
+            // Safety floors
+            if (_mass < 1e-9f)
+            {
+                _mass = 1e-9f;
+                _internalEnergy = _mass * Cv * 300f;
+                _airFraction = 1f;
+            }
+            if (_internalEnergy < 0f)
+                _internalEnergy = _mass * Cv * 300f;
+
+            // Final state update
+            if (_mass > 1e-12f && Volume > 1e-12f)
+            {
+                Temperature = _internalEnergy / (_mass * Cv);
+                Pressure = _mass * Rgas * Temperature / Volume;
+            }
+            else
+            {
+                Temperature = 300f;
+                Pressure = 101325f;
+            }
+
+            // Safety limits (unchanged, but now rarely triggered)
+            const float safetyPressure = 1.0f;
+            if (Pressure < safetyPressure && _mass > 1e-12f && Volume > 1e-12f)
+            {
+                Temperature = safetyPressure * Volume / (_mass * Rgas);
+                _internalEnergy = _mass * Cv * Temperature;
+                Pressure = safetyPressure;
+            }
+
+            const float maxPressure = 5e5f;
+            if (Pressure > maxPressure && _mass > 1e-12f && Volume > 1e-12f)
+            {
+                float targetMass = maxPressure * Volume / (Rgas * Temperature);
+                if (_mass > targetMass)
+                {
+                    _mass = targetMass;
+                    _internalEnergy = _mass * Cv * Temperature;
+                    Pressure = maxPressure;
+                }
+            }
+        }
+    }
+}

Components/Crankshaft.cs

@@ -4,33 +4,31 @@ namespace FluidSim.Components
 {
     public class Crankshaft
     {
-        public float AngularVelocity;   // rad/s
+        public float AngularVelocity;
-        public float CrankAngle;        // rad, 0 … 4π
+        public float CrankAngle;
         public float PreviousAngle;
 
-        public float Inertia = 0.2f;     // kg·m²
+        public float Inertia = 0.2f;
-        public float FrictionConstant;   // N·m
+        public float FrictionConstant;
-        public float FrictionViscous;    // N·m per rad/s
+        public float FrictionViscous;
-
         public float LastNetTorque { get; private set; }
-        public float AveragePower { get; private set; }     // smoothed, watts
+        public float AveragePower { get; private set; }
-        public float AverageTorque { get; private set; }    // smoothed, Nm
+        public float AverageTorque { get; private set; }
 
         private float externalTorque;
-        private float _loadTorque;       // external brake torque (Nm)
+        private float _loadTorque;
 
-        // Power averaging buffer
         private readonly float[] _powerBuffer;
-        private int _powerBufIdx;
+        private int _powerBufIdx, _powerBufCount;
-        private int _powerBufCount;
         private float _powerBufSum;
-
-        // Torque averaging buffer (same size as power buffer)
         private readonly float[] _torqueBuffer;
-        private int _torqueBufIdx;
+        private int _torqueBufIdx, _torqueBufCount;
-        private int _torqueBufCount;
         private float _torqueBufSum;
 
+        /// <summary>Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke.</summary>
+        public float CycleLength { get; set; } = 4f * MathF.PI;
+        public float CrankAngleRad => CrankAngle;
+
         public Crankshaft(float initialRPM = 400f)
         {
             AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
@@ -43,9 +41,13 @@ namespace FluidSim.Components
 
         public void AddTorque(float torque) => externalTorque += torque;
 
-        public void SetLoadTorque(float torque)
+        public void SetLoadTorque(float torque) => _loadTorque = Math.Max(torque, 0f);
+
+        private float _effectiveInertia;   // if >0, overrides Inertia
+
+        public void SetEffectiveInertia(float inertia)
         {
-            _loadTorque = Math.Max(torque, 0f);
+            _effectiveInertia = inertia;
         }
 
         public void Step(float dt)
@@ -57,51 +59,40 @@ namespace FluidSim.Components
 
             PreviousAngle = CrankAngle;
 
-            // Internal friction torque
             float friction = FrictionConstant * MathF.Sign(AngularVelocity)
                            + FrictionViscous * AngularVelocity;
 
-            // Net torque from gas pressure minus friction (used for power/torque display)
             float netTorque = externalTorque - friction;
             LastNetTorque = netTorque;
 
-            // Total torque after subtracting external load (brake)
             float totalNetTorque = netTorque - _loadTorque;
-            float alpha = totalNetTorque / Inertia;
+            float currentInertia = _effectiveInertia > 0f ? _effectiveInertia : Inertia;
+            float alpha = totalNetTorque / currentInertia;
             AngularVelocity += alpha * dt;
-
             if (AngularVelocity < 0f) AngularVelocity = 0f;
 
             CrankAngle += AngularVelocity * dt;
-            if (CrankAngle >= 4f * MathF.PI)
+            if (CrankAngle >= CycleLength)
-                CrankAngle -= 4f * MathF.PI;
+                CrankAngle -= CycleLength;
             else if (CrankAngle < 0f)
-                CrankAngle += 4f * MathF.PI;
+                CrankAngle += CycleLength;
 
-            // ---- Power averaging ----
+            // Power averaging
             float instantPower = netTorque * AngularVelocity;
             if (_powerBufCount == _powerBuffer.Length)
-            {
                 _powerBufSum -= _powerBuffer[_powerBufIdx];
-            }
             else
-            {
                 _powerBufCount++;
-            }
             _powerBuffer[_powerBufIdx] = instantPower;
             _powerBufSum += instantPower;
             _powerBufIdx = (_powerBufIdx + 1) % _powerBuffer.Length;
             AveragePower = _powerBufSum / _powerBufCount;
 
-            // ---- Torque averaging ----
+            // Torque averaging
             if (_torqueBufCount == _torqueBuffer.Length)
-            {
                 _torqueBufSum -= _torqueBuffer[_torqueBufIdx];
-            }
             else
-            {
                 _torqueBufCount++;
-            }
             _torqueBuffer[_torqueBufIdx] = netTorque;
             _torqueBufSum += netTorque;
             _torqueBufIdx = (_torqueBufIdx + 1) % _torqueBuffer.Length;

Components/Cylinder.cs

@@ -1,99 +1,25 @@
 using System;
-using System.Collections.Generic;
+using FluidSim.Components;   // if needed
-using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
-    public class Cylinder : IComponent
+    public class Cylinder : EngineCylinder
     {
-        public Port IntakePort { get; }
+        public float IVO, IVC, EVO, EVC;   // degrees in a 720° cycle
-        public Port ExhaustPort { get; }
-        public Crankshaft Crankshaft { get; }
 
-        private readonly Port[] _ports;
+        protected override float CycleLengthRad => 4f * MathF.PI;
-        IReadOnlyList<Port> IComponent.Ports => _ports;
+        protected override float MaxCycleDeg => 720f;
 
-        public float Bore { get; }
+        public override float IntakeValveArea =>
-        public float Stroke { get; }
+            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-        public float ConRodLength { get; }
+        public override float ExhaustValveArea =>
-        public float CompressionRatio { get; }
+            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-
-        public float IVO, IVC, EVO, EVC;     // degrees
-        public float IntakeValveDiameter = 0.03f;
-        public float ExhaustValveDiameter = 0.028f;
-        public float IntakeValveLift = 0.005f;
-        public float ExhaustValveLift = 0.005f;
-
-        public float IntakeValveMaxArea => MathF.PI * IntakeValveDiameter * IntakeValveLift;
-        public float ExhaustValveMaxArea => MathF.PI * ExhaustValveDiameter * ExhaustValveLift;
-
-        public float SparkAdvance = 20f;
-        public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f;
-        public float StoichiometricAFR = 14.7f;
-        public float FuelLowerHeatingValue = 44e6f;
-        public float EnergyVariationFraction = 0.05f;
-        public float MisfireProbability = 0.0f;
-        public float CylinderWallArea = 0.02f;
-        public float HeatTransferCoefficient = 100f;
-        public float AmbientTemperature = 300f;
-
-        public float PhaseOffset;   // rad
-
-        public float Volume => cylinderVolume;
-        public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
-        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
-        public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
-        public float Mass => _airMass + _exhaustMass;
-        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
-        public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
-
-        private float cylinderVolume, cylinderEnergy;
-        private float _airMass, _exhaustMass;
-        private float trappedAirMass, fuelMass, burnFraction;
-        private bool combustionActive, fuelInjected;
-        private float _energyFactor = 1f;
-        private readonly Random _random = new Random();
-
-        private const float Gamma = 1.4f;
-        private const float GasConstant = 287f;
-        private const float MaxPressurePa = 200e5f;
-        private const float MaxTemperatureK = 3500f;
 
         public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio,
                         float ivo, float ivc, float evo, float evc, Crankshaft crankshaft)
+            : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
         {
-            Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
-            CompressionRatio = compressionRatio;
             IVO = ivo; IVC = ivc; EVO = evo; EVC = evc;
-            Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
-
-            cylinderVolume = clearanceVolume;
-            float initRho = 1.225f;
-            _airMass = initRho * clearanceVolume;
-            _exhaustMass = 0f;
-            cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
-
-            IntakePort = new Port { Owner = this };
-            ExhaustPort = new Port { Owner = this };
-            _ports = new[] { IntakePort, ExhaustPort };
-        }
-
-        private float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
-        private float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
-        private float CrankRadius => Stroke * 0.5f;
-        private float Obliquity => CrankRadius / ConRodLength;
-
-        private float CrankDeg =>
-            ((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f;
-
-        public float ComputeVolume(float thetaRad)
-        {
-            float r = CrankRadius, l = ConRodLength;
-            float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
-            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
-            float x = r * (1f - cosTh) + l * (1f - term);
-            float area = MathF.PI * 0.25f * Bore * Bore;
-            return clearanceVolume + area * x;
         }
 
         private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
@@ -101,15 +27,10 @@ namespace FluidSim.Components
             float deg = thetaDeg % 720f;
             if (deg < 0f) deg += 720f;
 
-            float duration;
             float effectiveOpen = opens;
             float effectiveClose = closes;
-
+            if (closes < opens) effectiveClose += 720f;
-            if (closes < opens)
+            float duration = effectiveClose - effectiveOpen;
-            {
-                effectiveClose += 720f;
-            }
-            duration = effectiveClose - effectiveOpen;
             if (duration <= 0f) return 0f;
 
             float mapped = deg;
@@ -136,43 +57,9 @@ namespace FluidSim.Components
             return 0f;
         }
 
-        public float IntakeValveArea =>
+        protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
-            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-        public float ExhaustValveArea =>
-            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-
-        private float Wiebe(float angleSinceSpark)
         {
-            if (angleSinceSpark < WiebeStart) return 0f;
+            // Intake closing → fuel injection
-            float phi = (angleSinceSpark - WiebeStart) / WiebeDuration;
-            if (phi <= 0f) return 0f;
-            return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f));
-        }
-
-        public void PreStep(float dt)
-        {
-            // Speed‑dependent spark advance (simple linear)
-            float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-            SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);  // 10° at idle, ~30° at 10k rpm
-
-            float prevVolume = cylinderVolume;
-            float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
-            cylinderVolume = ComputeVolume(crankAngleRad);
-
-            float dV = cylinderVolume - prevVolume;
-            float pRel = Pressure - 101325f;
-            float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
-            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
-            float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
-            float pistonArea = MathF.PI * 0.25f * Bore * Bore;
-            Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
-
-            cylinderEnergy -= Pressure * dV;
-
-            float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % 720f;
-            float currDeg = crankAngleRad * 180f / MathF.PI % 720f;
-
-            // Intake closing – triggers fuel injection
             if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
             {
                 trappedAirMass = _airMass;
@@ -180,11 +67,14 @@ namespace FluidSim.Components
                 fuelInjected = true;
             }
 
-            // Spark
+            // Spark – occurs at TDC (0°) minus advance, every 720°
-            float sparkAngle = 0f - SparkAdvance;
+            float sparkAngle = (0f - SparkAdvance + 720f) % 720f;
-            if (sparkAngle < 0f) sparkAngle += 720f;
+            bool crossedSpark = false;
-            bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) ||
+            if (prevDeg < sparkAngle && currDeg >= sparkAngle)
-                                (prevDeg > sparkAngle + 360f && currDeg < sparkAngle);
+                crossedSpark = true;
+            else if (prevDeg > sparkAngle && currDeg < sparkAngle)
+                crossedSpark = true;
+
             if (crossedSpark && !combustionActive && fuelInjected)
             {
                 if (_random.NextDouble() < MisfireProbability)
@@ -199,7 +89,7 @@ namespace FluidSim.Components
                 }
             }
 
-            // Combustion
+            // Combustion progression
             if (combustionActive)
             {
                 float angleSinceSpark = currDeg - sparkAngle;
@@ -222,62 +112,6 @@ namespace FluidSim.Components
                     burnFraction = newFraction;
                 }
             }
-
-            // Heat loss
-            float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
-                            (Temperature - AmbientTemperature) * dt;
-            cylinderEnergy -= dQ_loss;
-
-            // Update port states
-            float p = Pressure, rho = Density, T = Temperature;
-            float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
-            float af = AirFraction;
-            IntakePort.Pressure = p; IntakePort.Density = rho;
-            IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
-            ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
-            ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
-        }
-
-        public void UpdateState(float dt)
-        {
-            float dmAir = 0f, dmExhaust = 0f, dE = 0f;
-            foreach (var port in _ports)
-            {
-                float mdot = port.MassFlowRate;
-                float af = mdot >= 0f ? port.AirFraction : AirFraction;
-                dmAir += mdot * af * dt;
-                dmExhaust += mdot * (1f - af) * dt;
-                dE += mdot * port.SpecificEnthalpy * dt;
-            }
-
-            _airMass += dmAir; _exhaustMass += dmExhaust;
-            cylinderEnergy += dE;
-
-            float V = MathF.Max(cylinderVolume, 1e-12f);
-            float currentP = (Gamma - 1f) * cylinderEnergy / V;
-            if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f);
-
-            float currentRho = (_airMass + _exhaustMass) / V;
-            float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f);
-            if (currentT > MaxTemperatureK)
-            {
-                float pAtTlimit = currentRho * GasConstant * MaxTemperatureK;
-                cylinderEnergy = pAtTlimit * V / (Gamma - 1f);
-            }
-
-            float totalMass = _airMass + _exhaustMass;
-            if (totalMass < 1e-9f)
-            {
-                _airMass = 1e-9f; _exhaustMass = 0f;
-                cylinderEnergy = 101325f * V / (Gamma - 1f);
-            }
-            else if (cylinderEnergy < 0f)
-            {
-                cylinderEnergy = 101325f * V / (Gamma - 1f);
-            }
-
-            if (_airMass < 0f) _airMass = 0f;
-            if (_exhaustMass < 0f) _exhaustMass = 0f;
         }
     }
 }

Components/EngineCylinder.cs

@@ -0,0 +1,203 @@
+using System;
+using System.Collections.Generic;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    /// <summary>Common base for all reciprocating engine cylinders.</summary>
+    public abstract class EngineCylinder : IComponent
+    {
+        public Port IntakePort { get; }
+        public Port ExhaustPort { get; }
+        public Crankshaft Crankshaft { get; }
+
+        private readonly Port[] _ports;
+        IReadOnlyList<Port> IComponent.Ports => _ports;
+
+        // ----- Geometry -----
+        public float Bore { get; }
+        public float Stroke { get; }
+        public float ConRodLength { get; }
+        public float CompressionRatio { get; }
+
+        // ----- Valve / port sizes (used for curtain area) -----
+        public float IntakeValveDiameter = 0.03f;
+        public float ExhaustValveDiameter = 0.028f;
+        public float IntakeValveLift = 0.005f;
+        public float ExhaustValveLift = 0.005f;
+
+        // ----- Combustion -----
+        public float SparkAdvance = 20f;
+        public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f;
+        public float StoichiometricAFR = 14.7f;
+        public float FuelLowerHeatingValue = 44e6f;
+        public float EnergyVariationFraction = 0.05f;
+        public float MisfireProbability = 0f;
+        public float CylinderWallArea = 0.02f;
+        public float HeatTransferCoefficient = 100f;
+        public float AmbientTemperature = 300f;
+
+        public float PhaseOffset;   // radians
+
+        // ----- State (public, used by drawing) -----
+        public float Volume => cylinderVolume;
+        public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
+        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
+        public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
+        public float Mass => _airMass + _exhaustMass;
+        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
+        public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
+
+        protected float cylinderVolume, cylinderEnergy;
+        protected float _airMass, _exhaustMass;
+        protected float trappedAirMass, fuelMass, burnFraction;
+        protected bool combustionActive, fuelInjected;
+        protected float _energyFactor = 1f;
+        protected readonly Random _random = new Random();
+
+        protected const float Gamma = 1.4f;
+        protected const float GasConstant = 287f;
+        protected const float MaxPressurePa = 200e5f;
+        protected const float MaxTemperatureK = 3500f;
+
+        // ----- Derived geometry (cycle‑independent) -----
+        protected float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
+        protected float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
+        protected float CrankRadius => Stroke * 0.5f;
+        protected float Obliquity => CrankRadius / ConRodLength;
+
+        // ----- Abstract members (cycle‑specific) -----
+        protected abstract float CycleLengthRad { get; }      // 4π or 2π
+        protected abstract float MaxCycleDeg { get; }         // 720 or 360
+        public abstract float IntakeValveArea { get; }
+        public abstract float ExhaustValveArea { get; }
+        protected abstract void HandleCycleEvents(float prevDeg, float currDeg, float dt);
+
+        protected EngineCylinder(float bore, float stroke, float conRodLength,
+                                 float compressionRatio, Crankshaft crankshaft)
+        {
+            Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
+            CompressionRatio = compressionRatio;
+            Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
+
+            cylinderVolume = clearanceVolume;
+            float initRho = 1.225f;
+            _airMass = initRho * clearanceVolume;
+            _exhaustMass = 0f;
+            cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
+
+            IntakePort = new Port { Owner = this };
+            ExhaustPort = new Port { Owner = this };
+            _ports = new[] { IntakePort, ExhaustPort };
+
+            // Set crankshaft cycle length
+            crankshaft.CycleLength = CycleLengthRad;
+        }
+
+        public float ComputeVolume(float thetaRad)
+        {
+            float r = CrankRadius, l = ConRodLength;
+            float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
+            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
+            float x = r * (1f - cosTh) + l * (1f - term);
+            float area = MathF.PI * 0.25f * Bore * Bore;
+            return clearanceVolume + area * x;
+        }
+
+        protected float CrankDeg =>
+            ((Crankshaft.CrankAngle + PhaseOffset) % CycleLengthRad) * 180f / MathF.PI;
+
+        protected float Wiebe(float angleSinceSpark)
+        {
+            if (angleSinceSpark < WiebeStart) return 0f;
+            float phi = (angleSinceSpark - WiebeStart) / WiebeDuration;
+            return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f));
+        }
+
+        // ----- Main update called before flow solver -----
+        public void PreStep(float dt)
+        {
+            // Speed‑dependent spark advance
+            float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);
+
+            float prevVolume = cylinderVolume;
+            float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
+            cylinderVolume = ComputeVolume(crankAngleRad);
+
+            // Piston work
+            float dV = cylinderVolume - prevVolume;
+            float pRel = Pressure - 101325f;
+            float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
+            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
+            float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
+            float pistonArea = MathF.PI * 0.25f * Bore * Bore;
+            Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
+
+            cylinderEnergy -= Pressure * dV;
+
+            float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % MaxCycleDeg;
+            float currDeg = crankAngleRad * 180f / MathF.PI % MaxCycleDeg;
+
+            // Let derived class handle valve events, spark, fuel
+            HandleCycleEvents(prevDeg, currDeg, dt);
+
+            // Heat loss
+            float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
+                            (Temperature - AmbientTemperature) * dt;
+            cylinderEnergy -= dQ_loss;
+
+            // Update port states
+            float p = Pressure, rho = Density, T = Temperature;
+            float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
+            float af = AirFraction;
+            IntakePort.Pressure = p; IntakePort.Density = rho;
+            IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
+            ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
+            ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
+        }
+
+        // ----- State update (mass/energy balance) -----
+        public void UpdateState(float dt)
+        {
+            float dmAir = 0f, dmExhaust = 0f, dE = 0f;
+            foreach (var port in _ports)
+            {
+                float mdot = port.MassFlowRate;
+                float af = mdot >= 0f ? port.AirFraction : AirFraction;
+                dmAir += mdot * af * dt;
+                dmExhaust += mdot * (1f - af) * dt;
+                dE += mdot * port.SpecificEnthalpy * dt;
+            }
+
+            _airMass += dmAir; _exhaustMass += dmExhaust;
+            cylinderEnergy += dE;
+
+            float V = MathF.Max(cylinderVolume, 1e-12f);
+            float currentP = (Gamma - 1f) * cylinderEnergy / V;
+            if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f);
+
+            float currentRho = (_airMass + _exhaustMass) / V;
+            float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f);
+            if (currentT > MaxTemperatureK)
+            {
+                float pAtTlimit = currentRho * GasConstant * MaxTemperatureK;
+                cylinderEnergy = pAtTlimit * V / (Gamma - 1f);
+            }
+
+            float totalMass = _airMass + _exhaustMass;
+            if (totalMass < 1e-9f)
+            {
+                _airMass = 1e-9f; _exhaustMass = 0f;
+                cylinderEnergy = 101325f * V / (Gamma - 1f);
+            }
+            else if (cylinderEnergy < 0f)
+            {
+                cylinderEnergy = 101325f * V / (Gamma - 1f);
+            }
+
+            if (_airMass < 0f) _airMass = 0f;
+            if (_exhaustMass < 0f) _exhaustMass = 0f;
+        }
+    }
+}

Components/TwoStrokeCylinder.cs

@@ -0,0 +1,192 @@
+// ============================================================
+// File: TwoStrokeCylinder.cs
+// ============================================================
+using System;
+using FluidSim.Interfaces;
+using FluidSim.Components;   // for Crankcase (if in same namespace)
+
+namespace FluidSim.Components
+{
+    /// <summary>
+    /// Two‑stroke cylinder with forced symmetrical port timings around BDC (180°).
+    /// Uses crankcase back‑pressure for accurate pumping work.
+    /// </summary>
+    public class TwoStrokeCylinder : EngineCylinder
+    {
+        // --- Port timing (computed from durations) ---
+        public float IVO => 180f - transferDuration / 2f;
+        public float IVC => 180f + transferDuration / 2f;
+        public float EVO => 180f - exhaustDuration / 2f;
+        public float EVC => 180f + exhaustDuration / 2f;
+
+        private readonly float transferDuration;   // degrees
+        private readonly float exhaustDuration;    // degrees
+
+        // --- Crankcase reference ---
+        private Crankcase? _crankcase;
+
+        protected override float CycleLengthRad => 2f * MathF.PI;
+        protected override float MaxCycleDeg => 360f;
+
+        public override float IntakeValveArea =>
+            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
+        public override float ExhaustValveArea =>
+            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
+
+        public TwoStrokeCylinder(float bore, float stroke, float conRodLength,
+                                 float compressionRatio,
+                                 float transferDuration, float exhaustDuration,
+                                 Crankshaft crankshaft)
+            : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
+        {
+            this.transferDuration = transferDuration;
+            this.exhaustDuration  = exhaustDuration;
+
+            if (EVO >= IVO)
+                throw new ArgumentException("Exhaust must open before transfer port.");
+        }
+
+        public void SetCrankcase(Crankcase crankcase)
+        {
+            _crankcase = crankcase;
+        }
+
+        // ----- Valve lift -----
+        private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
+        {
+            float deg = thetaDeg % 360f;
+            if (deg < 0f) deg += 360f;
+
+            float effectiveOpen = opens;
+            float effectiveClose = closes;
+            if (closes < opens) effectiveClose += 360f;
+            float duration = effectiveClose - effectiveOpen;
+            if (duration <= 0f) return 0f;
+
+            float mapped = deg;
+            if (mapped < opens) mapped += 360f;
+            if (mapped < opens || mapped > effectiveClose) return 0f;
+
+            float rampDur = duration * 0.25f;
+            float holdDur = duration - 2f * rampDur;
+
+            if (mapped >= opens && mapped < opens + rampDur)
+            {
+                float t = (mapped - opens) / rampDur;
+                return peakLift * t * t * (3f - 2f * t);
+            }
+            else if (mapped >= opens + rampDur && mapped < opens + rampDur + holdDur)
+            {
+                return peakLift;
+            }
+            else if (mapped >= opens + rampDur + holdDur && mapped <= effectiveClose)
+            {
+                float t = (mapped - (opens + rampDur + holdDur)) / rampDur;
+                return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
+            }
+            return 0f;
+        }
+
+        protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
+        {
+            // Transfer port closing → fuel injection
+            if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
+            {
+                trappedAirMass = _airMass;
+                fuelMass = trappedAirMass / StoichiometricAFR;
+                fuelInjected = true;
+            }
+
+            // Spark every 360° at TDC (0°) minus advance
+            float sparkAngle = (0f - SparkAdvance + 360f) % 360f;
+            bool crossedSpark = false;
+            if (prevDeg < sparkAngle && currDeg >= sparkAngle)
+                crossedSpark = true;
+            else if (prevDeg > sparkAngle && currDeg < sparkAngle)
+                crossedSpark = true;
+
+            if (crossedSpark && !combustionActive && fuelInjected)
+            {
+                if (_random.NextDouble() < MisfireProbability)
+                {
+                    combustionActive = false;
+                }
+                else
+                {
+                    combustionActive = true; burnFraction = 0f;
+                    float range = EnergyVariationFraction;
+                    _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
+                }
+            }
+
+            if (combustionActive)
+            {
+                float angleSinceSpark = currDeg - sparkAngle;
+                if (angleSinceSpark < 0f) angleSinceSpark += 360f;
+                float newFraction = Wiebe(angleSinceSpark);
+                if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance))
+                {
+                    newFraction = 1f; combustionActive = false;
+                    float totalMass = _airMass + _exhaustMass;
+                    _airMass = 0f; _exhaustMass = totalMass;
+                }
+                fuelInjected = false;
+
+                float dFraction = newFraction - burnFraction;
+                if (dFraction > 0f)
+                {
+                    float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction;
+                    cylinderEnergy += dQ;
+                    _exhaustMass += fuelMass * dFraction;
+                    burnFraction = newFraction;
+                }
+            }
+        }
+
+        // ----- Override torque calculation to use crankcase back‑pressure -----
+        public new void PreStep(float dt)
+        {
+            // Speed‑dependent spark advance
+            float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);
+
+            float prevVolume = cylinderVolume;
+            float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
+            cylinderVolume = ComputeVolume(crankAngleRad);
+
+            float dV = cylinderVolume - prevVolume;
+
+            // Use crankcase pressure as back‑pressure, ambient if not set
+            float backPressure = _crankcase?.Pressure ?? 101325f;
+            float pRel = Pressure - backPressure;
+
+            float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
+            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
+            float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
+            float pistonArea = MathF.PI * 0.25f * Bore * Bore;
+            Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
+
+            cylinderEnergy -= Pressure * dV;
+
+            float cycleLenDeg = 360f;
+            float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % cycleLenDeg;
+            float currDeg = crankAngleRad * 180f / MathF.PI % cycleLenDeg;
+
+            HandleCycleEvents(prevDeg, currDeg, dt);
+
+            // Heat loss
+            float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
+                            (Temperature - AmbientTemperature) * dt;
+            cylinderEnergy -= dQ_loss;
+
+            // Update port states
+            float p = Pressure, rho = Density, T = Temperature;
+            float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
+            float af = AirFraction;
+            IntakePort.Pressure = p; IntakePort.Density = rho;
+            IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
+            ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
+            ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
+        }
+    }
+}

Components/Vehicle.cs

@@ -0,0 +1,166 @@
+using System;
+
+namespace FluidSim.Components
+{
+    public class Vehicle
+    {
+        // ---- Gearbox ----
+        public int CurrentGear { get; private set; } = 0;
+        public readonly float[] GearRatios = { 2.5f, 1.8f, 1.4f, 1.1f, 0.9f, 0.75f };
+        public float FinalDriveRatio = 3.0f;
+        public float PrimaryReduction = 2.5f;
+
+        // ---- Clutch ----
+        public float ClutchInput { get; set; }
+        public float ClutchDisengageTime = 0.15f;
+        private float _clutchTimer;
+        private float _currentEngagement = 0f;
+
+        /// <summary>Time constant for clutch engagement smoothing (seconds).</summary>
+        public float EngagementSmoothTime = 0.5f;   // longer, gentler bite
+
+        private float TargetEngagement
+        {
+            get
+            {
+                if (ClutchInput > 0.01f) return 1f - ClutchInput;
+                if (CurrentGear == 0 || _clutchTimer > 0f) return 0f;
+                return 1f;
+            }
+        }
+
+        public float Engagement => _currentEngagement;
+
+        // ---- Clutch torque model ----
+        /// <summary>Peak clutch friction torque (Nm) when fully engaged at high RPM.</summary>
+        public float BaseMaxTorque = 80f;            // much lower than before
+
+        /// <summary>Stiffness when slipping (Nm per rad/s). Lower = softer engagement.</summary>
+        public float ClutchStiffness = 50f;          // very soft
+
+        /// <summary>Below this engine RPM, the clutch torque is progressively reduced to prevent stalling.</summary>
+        public float IdleRpm = 1200f;
+        public float StallPreventionRamp = 300f;     // RPM band above idle where torque ramps up
+
+        // ---- Physical constants ----
+        public float Mass = 160f;
+        public float WheelRadius = 0.32f;
+        public float DragCoefficient = 0.35f;
+        public float FrontalArea = 0.8f;
+        public float AirDensity = 1.225f;
+        public float RollingFrictionCoeff = 0.01f;
+        public float Gravity = 9.81f;
+
+        // ---- State ----
+        public float Speed { get; private set; }
+
+        public (float clutchTorqueOnEngine, float effectiveEngineInertia) Update(float engineRpm, float engineInertia, float dt)
+        {
+            if (_clutchTimer > 0f)
+            {
+                _clutchTimer -= dt;
+                if (_clutchTimer < 0f) _clutchTimer = 0f;
+            }
+
+            float target = TargetEngagement;
+            float smoothing = 1f - MathF.Exp(-dt / Math.Max(EngagementSmoothTime, 0.001f));
+            _currentEngagement += (target - _currentEngagement) * smoothing;
+            if (MathF.Abs(_currentEngagement - target) < 0.001f)
+                _currentEngagement = target;
+
+            float engagement = _currentEngagement;
+
+            float totalGear = 1f;
+            if (CurrentGear > 0)
+                totalGear = GearRatios[CurrentGear - 1] * FinalDriveRatio * PrimaryReduction;
+
+            float engineRadPerSec = engineRpm * 2f * MathF.PI / 60f;
+
+            float v = MathF.Max(Speed, 0f);
+            float drag = 0.5f * AirDensity * DragCoefficient * FrontalArea * v * v;
+            float rolling = RollingFrictionCoeff * Mass * Gravity;
+            float resistanceForce = drag + rolling;
+
+            float clutchTorque = 0f;
+            float effectiveInertia = engineInertia;
+
+            if (engagement > 0f && CurrentGear > 0)
+            {
+                float vehicleReflectedRadPerSec = (Speed / WheelRadius) * totalGear;
+                float slip = engineRadPerSec - vehicleReflectedRadPerSec;
+
+                // Stall prevention: reduce max torque when engine RPM is near idle
+                float torqueLimit = BaseMaxTorque * engagement;
+                if (engineRpm < IdleRpm + StallPreventionRamp)
+                {
+                    float factor = Math.Clamp((engineRpm - IdleRpm) / StallPreventionRamp, 0f, 1f);
+                    torqueLimit *= factor;
+                }
+
+                float stiffnessTorque = ClutchStiffness * engagement * slip;
+                clutchTorque = Math.Clamp(stiffnessTorque, -torqueLimit, torqueLimit);
+
+                // Lock if slip negligible and engagement high
+                if (engagement >= 0.99f && MathF.Abs(slip) < 1.0f)
+                {
+                    float vehicleInertia = Mass * WheelRadius * WheelRadius;
+                    float reflectedVehicleInertia = vehicleInertia / (totalGear * totalGear);
+                    effectiveInertia = engineInertia + reflectedVehicleInertia;
+
+                    Speed = engineRadPerSec * WheelRadius / totalGear;
+                    float loadTorque = resistanceForce * WheelRadius / totalGear;
+                    return (loadTorque, effectiveInertia);
+                }
+            }
+
+            float driveTorqueAtWheel = clutchTorque * totalGear;
+            float driveForce = driveTorqueAtWheel / WheelRadius;
+            float netForce = driveForce - resistanceForce;
+            float acceleration = netForce / Mass;
+            Speed += acceleration * dt;
+            if (Speed < 0f) Speed = 0f;
+
+            return (clutchTorque, engineInertia);
+        }
+
+        public void ShiftUp()
+        {
+            if (CurrentGear < GearRatios.Length)
+            {
+                CurrentGear++;
+                AutoDisengageClutch();
+            }
+        }
+
+        public void ShiftDown()
+        {
+            if (CurrentGear > 1)
+            {
+                CurrentGear--;
+                AutoDisengageClutch();
+            }
+        }
+
+        public void SetNeutral()
+        {
+            CurrentGear = 0;
+            _clutchTimer = 0f;
+        }
+
+        public void SetFirstGear()
+        {
+            if (CurrentGear == 0)
+            {
+                CurrentGear = 1;
+                AutoDisengageClutch();
+            }
+        }
+
+        private void AutoDisengageClutch()
+        {
+            _clutchTimer = ClutchDisengageTime;
+        }
+
+        public float SpeedKmh => Speed * 3.6f;
+    }
+}

Program.cs

@@ -48,13 +48,17 @@ public class Program
     private static float _loadTarget = 0.0f;          // 0‑1
     private static float _loadCurrent = 0.0f;
 
+    private static float _clutchTarget = 0f;
+    private static float _clutchCurrent = 0f;
+    private static bool _cKeyHeld = false;
+
     private const int TargetMaxFill = (int)(SampleRate * 0.2);
 
     public static void Main()
     {
         var window = CreateWindow();
         LoadFont();
-        _scenario = new SingleCylScenario();
+        _scenario = new TwoStrokeScenario();
         _scenario.Font = _overlayFont;
         _scenario.Initialize(SampleRate);
         _lastThrottleUpdateTime = 0.0f;
@@ -102,6 +106,11 @@ public class Program
 
             _scenario.Throttle = _throttleCurrent;
 
+            float clutchDesired = _cKeyHeld ? 1f : 0f;
+            float clutchSmoothing = 1f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
+            _clutchCurrent += (clutchDesired - _clutchCurrent) * clutchSmoothing;
+            _scenario.Clutch = _clutchCurrent;
+
 
             // ---- Drawing ----
             if (now - lastDrawTime >= 1.0 / DrawFrequency)
@@ -111,6 +120,7 @@ public class Program
                     string toggleHint = _isRealTime ? "[Space] slow mo" : "[Space] real time";
                     _overlayText.DisplayedString =
                         $"{toggleHint}  Speed: {_currentDisplaySpeed:F3}x  RT: {(_currentDisplaySpeed * 100.0):F1}%  Sim load: {_loadTracker.LoadPercent:F0}%\n" +
+                        $"Clutch: {_clutchCurrent*100:F0}%  [C]" +
                         $"Load: {_loadCurrent*100:F0}%  [←][→]  Throttle: {_throttleCurrent * 100:F0}%  Target: {_throttleTarget * 100:F0}%  [W] {(_wKeyHeld ? "BLIP" : "---")}";
                 }
 
@@ -221,6 +231,17 @@ public class Program
             case Keyboard.Key.Right:
                 _loadTarget = MathF.Min(1.0f, _loadTarget + 0.05f);
                 break;
+
+            case Keyboard.Key.E:
+                _scenario.ShiftUp();
+                break;
+            case Keyboard.Key.Q:
+                _scenario.ShiftDown();
+                break;
+
+            case Keyboard.Key.C:
+                _cKeyHeld = true;
+                break;
         }
     }
 
@@ -228,5 +249,8 @@ public class Program
     {
         if (e.Code == Keyboard.Key.W)
             _wKeyHeld = false;
+
+        if (e.Code == Keyboard.Key.C)
+            _cKeyHeld = false;
     }
 }

Scenarios/Scenario.cs

@@ -13,12 +13,16 @@ namespace FluidSim.Tests
         protected const float AmbientTemperature = 300f;
         public float Throttle { get; set; }
         public float Load { get; set; }
+        public float Clutch { get; set; }   // 0 = engaged, 1 = fully disengaged (manual lever)
         public Font? Font { get; set; }
 
         public abstract void Initialize(int sampleRate);
         public abstract float Process();
         public abstract void Draw(RenderWindow target);
 
+        public virtual void ShiftUp() { }
+        public virtual void ShiftDown() { }     
+
         // ---- Dyno curve graph ----
         private const float RpmBinSize = 50f;
         private readonly List<(float powerKw, float torqueNm)> _dynoBins = new();
@@ -259,7 +263,7 @@ namespace FluidSim.Tests
             target.Draw(border);
         }
 
-        protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
+        protected void DrawCylinder(RenderWindow target, EngineCylinder cylinder,
                                     float centerX, float topY, float width, float maxHeight)
         {
             float fraction = cylinder.PistonFraction;
@@ -298,7 +302,8 @@ namespace FluidSim.Tests
         }
 
         protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex,
-                                float pipeCenterY, float pipeStartX, float pipeEndX)
+                        float pipeCenterY, float pipeStartX, float pipeEndX,
+                        float areaScale = 0f)
         {
             int start = pipeSystem.GetPipeStart(pipeIndex);
             int end = pipeSystem.GetPipeEnd(pipeIndex);
@@ -307,20 +312,34 @@ namespace FluidSim.Tests
 
             float pipeLen = pipeEndX - pipeStartX;
             float dx = pipeLen / (n - 1);
-            float baseRadius = 25f;
 
             var centers = new float[n];
             var radii = new float[n];
             var temps = new float[n];
+
             for (int i = 0; i < n; i++)
             {
                 int cell = start + i;
                 float p = pipeSystem.GetCellPressure(cell);
                 float rho = pipeSystem.GetCellDensity(cell);
                 temps[i] = p / MathF.Max(rho * 287f, 1e-12f);
+
+                if (areaScale > 0f)
+                {
+                    // Use actual cell area to determine visual radius
+                    float area = pipeSystem.GetCellArea(cell);
+                    radii[i] = MathF.Sqrt(area / MathF.PI) * areaScale;
+                    if (radii[i] < 1f) radii[i] = 1f;
+                }
+                else
+                {
+                    // Original pressure‑based radius
                     float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
+                    float baseRadius = 25f;   // default visual radius for constant‑area pipes
                     radii[i] = baseRadius * (1f + dev * 2f);
                     if (radii[i] < 2f) radii[i] = 2f;
+                }
+
                 centers[i] = pipeStartX + i * dx;
             }
 

Scenarios/TwoStrokeScenario.cs

@@ -0,0 +1,288 @@
+using FluidSim.Components;
+using FluidSim.Core;
+using FluidSim.Interfaces;
+using FluidSim.Utils;
+using SFML.Graphics;
+using SFML.System;
+using System;
+
+namespace FluidSim.Tests
+{
+    public class TwoStrokeScenario : Scenario
+    {
+        private Crankshaft crankshaft;
+        private TwoStrokeCylinder cylinder;
+        private Crankcase crankcase;
+
+        private PipeSystem pipeSystem;
+        private BoundarySystem boundaries;
+        private Solver solver;
+
+        private Volume0D intakePlenum;
+        private Port plenumInlet, plenumOutlet;
+
+        private int throttleAreaIdx, reedInletIdx, reedOutletIdx,
+                    transferInletIdx, transferOutletIdx, exhaustValveIdx;
+        private float[] orificeAreas;
+        private int intakeOpenIdx, exhaustOpenIdx;
+
+        private SoundProcessor exhaustSound, intakeSound;
+        private OutdoorExhaustReverb reverb;
+
+        private double dt;
+        private int stepCount;
+
+        private float maxThrottleArea;
+        private float intakePipeArea, reedPipeArea, transferPipeArea, exhaustHeaderArea;
+        private bool reedOpen;
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm);
+
+            // ---- Crankshaft ----
+            crankshaft = new Crankshaft(3000);
+            crankshaft.CycleLength = 2f * MathF.PI;
+            crankshaft.Inertia = 0.01f;
+            crankshaft.FrictionConstant = 1.0f;
+            crankshaft.FrictionViscous = 0.002f;
+
+            // ---- Cylinder (125cc) ----
+            float bore = 0.054f, stroke = 0.0545f, conRod = 0.110f, compRatio = 7.2f;
+            float transferDur = 140f, exhaustDur = 195f;
+            cylinder = new TwoStrokeCylinder(bore, stroke, conRod, compRatio,
+                                            transferDur, exhaustDur, crankshaft)
+            {
+                // FIX: realistic transfer port diameter (was 40mm)
+                IntakeValveDiameter  = 0.030f,   // 30 mm
+                IntakeValveLift      = 0.010f,
+                ExhaustValveDiameter = 0.040f,
+                ExhaustValveLift     = 0.010f
+            };
+
+            // ---- Crankcase ----
+            float crankRadius = stroke * 0.5f;
+            float ccClearance = 150e-6f;
+            crankcase = new Crankcase(crankshaft, crankRadius, conRod, bore,
+                                    ccClearance, 101325f, 300f);
+            cylinder.SetCrankcase(crankcase);
+
+            // ---- Pipe system ----
+            int intakeCells    = 8;
+            int reedCells      = 4;
+            int transferCells  = 8;
+            int exhaustCells   = 60;
+            int totalCells = intakeCells + reedCells + transferCells + exhaustCells;
+
+            int[] pipeStart = {
+                0,
+                intakeCells,
+                intakeCells + reedCells,
+                intakeCells + reedCells + transferCells
+            };
+            int[] pipeEnd = {
+                intakeCells,
+                intakeCells + reedCells,
+                intakeCells + reedCells + transferCells,
+                totalCells
+            };
+
+            float[] area = new float[totalCells];
+            float[] dx   = new float[totalCells];
+
+            float intakeDia = 0.042f, reedDia = 0.040f, transferDia = 0.040f;
+            intakePipeArea   = MathF.PI * 0.25f * intakeDia * intakeDia;
+            reedPipeArea     = MathF.PI * 0.25f * reedDia * reedDia;
+            transferPipeArea = MathF.PI * 0.25f * transferDia * transferDia;
+
+            for (int i = 0; i < intakeCells; i++)
+            { area[i] = intakePipeArea; dx[i] = 0.100f / intakeCells; }
+
+            for (int i = intakeCells; i < intakeCells + reedCells; i++)
+            { area[i] = reedPipeArea; dx[i] = 0.030f / reedCells; }
+
+            for (int i = intakeCells + reedCells; i < intakeCells + reedCells + transferCells; i++)
+            { area[i] = transferPipeArea; dx[i] = 0.200f / transferCells; }
+
+            float hdrD = 0.040f, hdrL = 0.130f;
+            float difEndD = 0.070f, difL = 0.250f;
+            float belL = 0.220f;
+            float convEndD = 0.028f, convL = 0.160f;
+            float stiL = 0.080f;
+            float totL = hdrL + difL + belL + convL + stiL;
+            exhaustHeaderArea = MathF.PI * 0.25f * hdrD * hdrD;
+            float bellyArea = MathF.PI * 0.25f * difEndD * difEndD;
+            float stingerArea = MathF.PI * 0.25f * convEndD * convEndD;
+
+            int exhStart = intakeCells + reedCells + transferCells;
+            int hdrC = (int)(exhaustCells * hdrL / totL);
+            int difC = (int)(exhaustCells * difL / totL);
+            int belC = (int)(exhaustCells * belL / totL);
+            int conC = (int)(exhaustCells * convL / totL);
+            int stiC = exhaustCells - hdrC - difC - belC - conC;
+
+            int idx = 0;
+            for (int i = exhStart; i < totalCells; i++)
+            {
+                if (idx < hdrC) { area[i] = exhaustHeaderArea; dx[i] = hdrL / hdrC; }
+                else if (idx < hdrC + difC)
+                {
+                    float t = (idx - hdrC) / (float)(difC - 1);
+                    float dia = hdrD + (difEndD - hdrD) * t;
+                    area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = difL / difC;
+                }
+                else if (idx < hdrC + difC + belC) { area[i] = bellyArea; dx[i] = belL / belC; }
+                else if (idx < hdrC + difC + belC + conC)
+                {
+                    float t = (idx - hdrC - difC - belC) / (float)(conC - 1);
+                    float dia = difEndD + (convEndD - difEndD) * t;
+                    area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = convL / conC;
+                }
+                else { area[i] = stingerArea; dx[i] = stiL / stiC; }
+                idx++;
+            }
+
+            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
+                                        1.225f, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 0.8f;
+            pipeSystem.EnergyRelaxationRate = 0.4f;
+
+            // ---- Volumes ----
+            intakePlenum = new Volume0D(0.5e-3f, 101325f, 300f);
+            plenumInlet  = intakePlenum.CreatePort();
+            plenumOutlet = intakePlenum.CreatePort();
+
+            // ---- Boundary system ----
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 6, maxOpenEnds: 2);
+            throttleAreaIdx   = 0;
+            reedInletIdx      = 1;
+            reedOutletIdx     = 2;
+            transferInletIdx  = 3;
+            transferOutletIdx = 4;
+            exhaustValveIdx   = 5;
+
+            boundaries.AddOpenEnd(0, true,  101325f, intakePipeArea);
+            intakeOpenIdx = 0;
+            boundaries.AddOpenEnd(3, false, 101325f, stingerArea);
+            exhaustOpenIdx = 1;
+
+            boundaries.AddOrifice(plenumInlet,           0, false, throttleAreaIdx, 0.72f);
+            boundaries.AddOrifice(plenumOutlet,          1, true,  reedInletIdx,   1.0f);
+            boundaries.AddOrifice(crankcase.IntakePort,  1, false, reedOutletIdx,  0.9f);
+            boundaries.AddOrifice(crankcase.TransferPort,2, true,  transferInletIdx,1.0f);
+            boundaries.AddOrifice(cylinder.IntakePort,   2, false, transferOutletIdx,1.0f);
+            boundaries.AddOrifice(cylinder.ExhaustPort,  3, true,  exhaustValveIdx, 0.7f);
+
+            orificeAreas = new float[6];
+            orificeAreas[reedInletIdx]      = reedPipeArea;
+            orificeAreas[reedOutletIdx]     = 0f;
+            orificeAreas[transferInletIdx]  = transferPipeArea;
+            orificeAreas[transferOutletIdx] = 0f;
+
+            // ---- Solver ----
+            solver = new Solver { SubStepCount = 4 };
+            solver.SetTimeStep(dt);
+            solver.SetPipeSystem(pipeSystem);
+            solver.SetBoundarySystem(boundaries);
+            solver.AddComponent(cylinder);
+            solver.AddComponent(crankcase);
+            solver.AddComponent(intakePlenum);
+
+            // ---- Sound ----
+            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
+            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
+            reverb = new OutdoorExhaustReverb(sampleRate);
+
+            stepCount = 0;
+            Console.WriteLine("Two‑Stroke engine ready.");
+        }
+
+        public override float Process()
+        {
+            const float reedMargin = 200f;
+            if (crankcase.Pressure < intakePlenum.Pressure - reedMargin)
+                reedOpen = true;
+            else if (crankcase.Pressure > intakePlenum.Pressure + reedMargin)
+                reedOpen = false;
+
+            float throttledFraction = Throttle;
+            if (throttledFraction < 0.001f) throttledFraction = 0f;
+            throttledFraction = Math.Clamp(throttledFraction, 0f, 1f);
+            float throttledArea = maxThrottleArea * throttledFraction;
+
+            orificeAreas[throttleAreaIdx] = throttledArea;
+            orificeAreas[reedOutletIdx] = reedOpen ? reedPipeArea : 0f;
+            orificeAreas[transferOutletIdx] = cylinder.IntakeValveArea;
+            orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
+            boundaries.SetOrificeAreas(orificeAreas);
+
+            if (stepCount < 20000)
+                crankshaft.AddTorque(5.0f);
+
+            // FIX: update crankshaft BEFORE volumes, so crankcase and cylinder see the same new angle
+            crankshaft.Step((float)dt);
+            cylinder.PreStep((float)dt);
+            crankcase.PreStep((float)dt);
+
+            solver.Step();
+            stepCount++;
+
+            float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
+            float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
+
+            float exhaustDry = exhaustSound.Process(exhaustFlow);
+            float intakeDry  = intakeSound.Process(intakeFlow);
+
+            if (stepCount % 2000 == 0)
+            {
+                float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+                Console.WriteLine($"Step {stepCount} | RPM={rpm:F0} | CylP={cylinder.Pressure/1e5f:F2} bar | CCP={crankcase.Pressure/1e5f:F3} bar | Plenum={intakePlenum.Pressure/1e5f:F3} bar | Reed={reedOpen}");
+            }
+
+            return reverb.Process((intakeDry + exhaustDry) * 0.5f);
+        }
+
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+
+            float startX = 40f;
+            float endX = winW - 80f;
+
+            DrawPipe(target, pipeSystem, 0, winH * 0.25f, startX, startX + 120f);
+            var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
+            {
+                FillColor = Color.Yellow,
+                Position = new Vector2f(startX + 125f, winH * 0.25f - 15f)
+            };
+            target.Draw(throttleRect);
+            float plenX = startX + 140f;
+            DrawVolume(target, intakePlenum, plenX + 30f, winH * 0.25f - 25f, 60f, 50f);
+
+            float reedStartX = plenX + 70f;
+            DrawPipe(target, pipeSystem, 1, winH * 0.25f, reedStartX, reedStartX + 30f);
+
+            float transStartX = reedStartX + 40f;
+            DrawPipe(target, pipeSystem, 2, winH * 0.45f, transStartX, transStartX + 120f);
+
+            float cylCX = transStartX + 180f;
+            float cylTopY = winH * 0.45f - 90f;
+            DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f);
+
+            float exhStartX = cylCX + 60f;
+            DrawPipe(target, pipeSystem, 3, winH * 0.65f, exhStartX, endX, areaScale: 800f);
+
+            float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            float powerKw = crankshaft.AveragePower * 1e-3f;
+            DrawLabel(target, $"RPM: {rpm:F0}", new Vector2f(20, 90), Color.White, 24);
+            DrawLabel(target, $"Power: {powerKw:F2} kW", new Vector2f(20, 115), Color.White, 24);
+
+            float torqueNm = crankshaft.AverageTorque;
+            UpdateDynoCurve(rpm, powerKw, torqueNm);
+            DrawDynoCurve(target, winW - 410f, winH - 260f, 400f, 250f, rpm, powerKw);
+        }
+    }
+}