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

@@ -27,6 +27,7 @@ namespace FluidSim.Components
 
         /// <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)
         {

Components/TwoStrokeCylinder.cs

@@ -1,124 +1,111 @@
+// ============================================================
+// 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 symmetrical port timings centred on BDC (180°).
+    /// Two‑stroke cylinder with forced symmetrical port timings around BDC (180°).
-    ///
+    /// Uses crankcase back‑pressure for accurate pumping work.
-    /// Changes vs. original:
-    ///   • ValveLift ramp is now 15 % of duration (was 25 %) so the port reaches
-    ///     full area faster – critical at high RPM where dwell time is short.
-    ///   • Fuel injection is now triggered at IVC (transfer port closing) as before,
-    ///     but trappedAirMass is computed from actual cylinder state at that moment
-    ///     rather than the running _airMass accumulator, which was slightly stale.
-    ///   • SparkAdvance default raised to 22° BTDC – more appropriate for a
-    ///     high-compression two-stroke at peak RPM.  The scenario can still override it.
     /// </summary>
     public class TwoStrokeCylinder : EngineCylinder
     {
-        // ── Port timing read-outs (degrees, 0 = TDC) ───────────────────────────
+        // --- Port timing (computed from durations) ---
-        public float IVO => 180f - TransferDuration / 2f;   // transfer opens
+        public float IVO => 180f - transferDuration / 2f;
-        public float IVC => 180f + TransferDuration / 2f;   // transfer closes
+        public float IVC => 180f + transferDuration / 2f;
-        public float EVO => 180f - ExhaustDuration  / 2f;   // exhaust opens
+        public float EVO => 180f - exhaustDuration / 2f;
-        public float EVC => 180f + ExhaustDuration  / 2f;   // exhaust closes
+        public float EVC => 180f + exhaustDuration / 2f;
 
-        // ── Configurable durations ──────────────────────────────────────────────
+        private readonly float transferDuration;   // degrees
-        public float TransferDuration { get; }   // default: 155°
+        private readonly float exhaustDuration;    // degrees
-        public float ExhaustDuration  { get; }   // default: 195°
 
-        // Fraction of port-open duration used for ramp-up / ramp-down.
+        // --- Crankcase reference ---
-        // 0.15 → port at full area for the middle 70 % of open time.
+        private Crankcase? _crankcase;
-        private const float RampFraction = 0.15f;
 
         protected override float CycleLengthRad => 2f * MathF.PI;
-        protected override float MaxCycleDeg    => 360f;
+        protected override float MaxCycleDeg => 360f;
 
         public override float IntakeValveArea =>
-            MathF.PI * IntakeValveDiameter
+            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-                     * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-
         public override float ExhaustValveArea =>
-            MathF.PI * ExhaustValveDiameter
+            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-                     * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
 
-        // ── Constructor ─────────────────────────────────────────────────────────
         public TwoStrokeCylinder(float bore, float stroke, float conRodLength,
                                  float compressionRatio,
                                  float transferDuration, float exhaustDuration,
                                  Crankshaft crankshaft)
             : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
         {
-            TransferDuration = transferDuration;
+            this.transferDuration = transferDuration;
-            ExhaustDuration  = exhaustDuration;
+            this.exhaustDuration  = exhaustDuration;
 
             if (EVO >= IVO)
-                throw new ArgumentException(
+                throw new ArgumentException("Exhaust must open before transfer port.");
-                    $"Exhaust must open before transfer port. " +
-                    $"EVO={EVO:F1}° must be less than IVO={IVO:F1}°. " +
-                    $"Increase exhaustDuration or decrease transferDuration.");
         }
 
-        // ── Valve lift profile ──────────────────────────────────────────────────
+        public void SetCrankcase(Crankcase crankcase)
-        /// <summary>
+        {
-        /// Smooth trapezoidal lift: fast ramp (15 % of duration), flat top (70 %),
+            _crankcase = crankcase;
-        /// fast ramp-down (15 %).  Ramps use a smoothstep (3t²-2t³) curve so the
+        }
-        /// area derivative is C1-continuous (no kink at ramp/plateau boundaries).
+
-        /// </summary>
+        // ----- Valve lift -----
-        private static float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
+        private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
         {
-            // Normalise to [0, 360)
             float deg = thetaDeg % 360f;
             if (deg < 0f) deg += 360f;
 
-            // Handle wrap-around (e.g. opens=170°, closes=190° is fine;
+            float effectiveOpen = opens;
-            // a port that crosses 360° would need closes+360).
+            float effectiveClose = closes;
-            float effectiveClose = closes < opens ? closes + 360f : closes;
+            if (closes < opens) effectiveClose += 360f;
-            float duration = effectiveClose - opens;
+            float duration = effectiveClose - effectiveOpen;
             if (duration <= 0f) return 0f;
 
-            // Map deg into the same number-line as opens/effectiveClose
+            float mapped = deg;
-            float mapped = deg < opens ? deg + 360f : deg;
+            if (mapped < opens) mapped += 360f;
             if (mapped < opens || mapped > effectiveClose) return 0f;
 
-            float rampDur = duration * RampFraction;
+            float rampDur = duration * 0.25f;
-            float holdEnd = effectiveClose - rampDur;
+            float holdDur = duration - 2f * rampDur;
 
-            if (mapped < opens + rampDur)
+            if (mapped >= opens && mapped < opens + rampDur)
             {
-                // Opening ramp: smoothstep
                 float t = (mapped - opens) / rampDur;
                 return peakLift * t * t * (3f - 2f * t);
             }
-            else if (mapped <= holdEnd)
+            else if (mapped >= opens + rampDur && mapped < opens + rampDur + holdDur)
             {
-                // Flat top – full area
                 return peakLift;
             }
-            else
+            else if (mapped >= opens + rampDur + holdDur && mapped <= effectiveClose)
             {
-                // Closing ramp: smoothstep reversed
+                float t = (mapped - (opens + rampDur + holdDur)) / rampDur;
-                float t = (mapped - holdEnd) / rampDur;
                 return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
             }
+            return 0f;
         }
 
-        // ── Cycle event handler ─────────────────────────────────────────────────
         protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
         {
-            // ── Fuel injection at transfer-port closing (IVC) ──────────────────
+            // Transfer port closing → fuel injection
-            // At IVC the cylinder is sealed; whatever air is trapped is what we burn.
+            if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
-            if (CrossedAngle(prevDeg, currDeg, IVC))
             {
                 trappedAirMass = _airMass;
-                fuelMass       = trappedAirMass / StoichiometricAFR;
+                fuelMass = trappedAirMass / StoichiometricAFR;
-                fuelInjected   = true;
+                fuelInjected = true;
             }
 
-            // ── Ignition ───────────────────────────────────────────────────────
+            // Spark every 360° at TDC (0°) minus advance
-            // SparkAdvance default is ~22° BTDC on the base class; scenario can override.
+            float sparkAngle = (0f - SparkAdvance + 360f) % 360f;
-            float sparkAngle = (360f - SparkAdvance) % 360f;
+            bool crossedSpark = false;
+            if (prevDeg < sparkAngle && currDeg >= sparkAngle)
+                crossedSpark = true;
+            else if (prevDeg > sparkAngle && currDeg < sparkAngle)
+                crossedSpark = true;
 
-            if (CrossedAngle(prevDeg, currDeg, sparkAngle) && !combustionActive && fuelInjected)
+            if (crossedSpark && !combustionActive && fuelInjected)
             {
                 if (_random.NextDouble() < MisfireProbability)
                 {
@@ -126,58 +113,80 @@ namespace FluidSim.Components
                 }
                 else
                 {
-                    combustionActive = true;
+                    combustionActive = true; burnFraction = 0f;
-                    burnFraction     = 0f;
                     float range = EnergyVariationFraction;
                     _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
                 }
             }
 
-            // ── Combustion heat release (Wiebe) ────────────────────────────────
             if (combustionActive)
             {
                 float angleSinceSpark = currDeg - sparkAngle;
                 if (angleSinceSpark < 0f) angleSinceSpark += 360f;
-
                 float newFraction = Wiebe(angleSinceSpark);
-                bool burnComplete = newFraction >= 1f
+                if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance))
-                                 || angleSinceSpark > WiebeDuration + WiebeStart + SparkAdvance;
-
-                if (burnComplete)
                 {
-                    newFraction      = 1f;
+                    newFraction = 1f; combustionActive = false;
-                    combustionActive = false;
+                    float totalMass = _airMass + _exhaustMass;
-                    fuelInjected     = false;
+                    _airMass = 0f; _exhaustMass = totalMass;
-                    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;
+                    _exhaustMass += fuelMass * dFraction;
-                    burnFraction    = newFraction;
+                    burnFraction = newFraction;
                 }
             }
         }
 
-        // ── Helper: did the crank cross a target angle this step? ───────────────
+        // ----- Override torque calculation to use crankcase back‑pressure -----
-        /// <summary>
+        public new void PreStep(float dt)
-        /// Returns true if the crank swept through <paramref name="target"/> going
-        /// from <paramref name="prev"/> to <paramref name="curr"/> in a single step.
-        /// Handles wrap-around at 360°.
-        /// </summary>
-        private static bool CrossedAngle(float prev, float curr, float target)
         {
-            // Normal case (no wrap)
+            // Speed‑dependent spark advance
-            if (curr >= prev)
+            float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-                return prev < target && target <= curr;
+            SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);
 
-            // Wrapped past 360° → two intervals to check
+            float prevVolume = cylinderVolume;
-            return prev < target || target <= curr;
+            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;
         }
     }
 }

Scenarios/TwoStrokeScenario.cs

@@ -12,6 +12,7 @@ namespace FluidSim.Tests
     {
         private Crankshaft crankshaft;
         private TwoStrokeCylinder cylinder;
+        private Crankcase crankcase;
 
         private PipeSystem pipeSystem;
         private BoundarySystem boundaries;
@@ -19,12 +20,9 @@ namespace FluidSim.Tests
 
         private Volume0D intakePlenum;
         private Port plenumInlet, plenumOutlet;
-        private Volume0D exhaustMuffler;
-        private Port mufflerIn, mufflerOut;
 
-        private Vehicle vehicle;
+        private int throttleAreaIdx, reedInletIdx, reedOutletIdx,
-
+                    transferInletIdx, transferOutletIdx, exhaustValveIdx;
-        private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
         private float[] orificeAreas;
         private int intakeOpenIdx, exhaustOpenIdx;
 
@@ -34,233 +32,199 @@ namespace FluidSim.Tests
         private double dt;
         private int stepCount;
 
-        private float _maxThrottleArea;
+        private float maxThrottleArea;
-        private float intakePipeArea, exhaustHeaderArea;
+        private float intakePipeArea, reedPipeArea, transferPipeArea, exhaustHeaderArea;
-
+        private bool reedOpen;
-        public override void ShiftUp()   => vehicle.ShiftUp();
-        public override void ShiftDown() => vehicle.ShiftDown();
 
         public override void Initialize(int sampleRate)
         {
             dt = 1.0 / sampleRate;
 
-            // ── Vehicle ──────────────────────────────────────────────────────────
+            maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm);
-            vehicle = new Vehicle();
 
-            // ── Throttle body: 42 mm – wider to reduce high-RPM intake restriction ──
+            // ---- Crankshaft ----
-            _maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm);
+            crankshaft = new Crankshaft(3000);
-
+            crankshaft.CycleLength = 2f * MathF.PI;
-            // ── Crankshaft ───────────────────────────────────────────────────────
+            crankshaft.Inertia = 0.01f;
-            // Lighter flywheel for quicker revving; friction tuned to ~0.5 kW loss at idle
+            crankshaft.FrictionConstant = 1.0f;
-            crankshaft = new Crankshaft(2000);
+            crankshaft.FrictionViscous = 0.002f;
-            crankshaft.CycleLength      = 2f * MathF.PI;   // two-stroke: fire every rev
-            crankshaft.Inertia          = 0.06f;            // lighter flywheel
-            crankshaft.FrictionConstant = 0.4f;             // ~0.4 Nm constant drag
-            crankshaft.FrictionViscous  = 0.0004f;          // ~2.5 Nm at 10 000 RPM
-
-            // ── Cylinder: 125 cc, motocross-style two-stroke ─────────────────────
-            // Bore × stroke = 54 × 54.5 mm → 124.9 cc
-            float bore         = 0.054f;
-            float stroke       = 0.0545f;
-            float conRod       = 0.110f;   // ~2× stroke
-            float compRatio    = 7.2f;     // geometric CR; effective CR after port closure is ~12:1
-
-            // Port timings: exhaust 195°, transfer 155° – competitive MX 125
-            float transferDuration = 155f;
-            float exhaustDuration  = 195f;
 
+            // ---- 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,
-                                             transferDuration, exhaustDuration,
+                                            transferDur, exhaustDur, crankshaft)
-                                             crankshaft)
             {
-                IntakeValveDiameter  = 0.042f,   // matched to intake pipe
+                // FIX: realistic transfer port diameter (was 40mm)
-                IntakeValveLift      = 0.015f,
+                IntakeValveDiameter  = 0.030f,   // 30 mm
+                IntakeValveLift      = 0.010f,
                 ExhaustValveDiameter = 0.040f,
-                ExhaustValveLift     = 0.013f
+                ExhaustValveLift     = 0.010f
             };
 
-            // ── Pipe geometry ────────────────────────────────────────────────────
+            // ---- Crankcase ----
-            //
+            float crankRadius = stroke * 0.5f;
-            // Layout (all lengths in mm):
+            float ccClearance = 150e-6f;
-            //   Intake path:  airbox stub 100 mm  |  runner 180 mm
+            crankcase = new Crankcase(crankshaft, crankRadius, conRod, bore,
-            //   Exhaust path: expansion chamber tuned to ~9 000 RPM power peak
+                                    ccClearance, 101325f, 300f);
-            //     header      170 mm  Ø 40 mm
+            cylinder.SetCrankcase(crankcase);
-            //     diffuser    280 mm  Ø 40 → 72 mm
-            //     belly       200 mm  Ø 72 mm
-            //     convergent  130 mm  Ø 72 → 28 mm
-            //     stinger      70 mm  Ø 28 mm
-            //     total       850 mm
-            //
-            // Cell sizing: ~14 mm/cell.
-            // CFL: c_sound ≈ 550 m/s, dx=0.014 m → dt_max ≈ 25 µs
-            //      at 44100 Hz dt = 22.7 µs → SubStepCount=4 keeps CFL safely ≤ 1
 
-            // --- Cell counts ---
+            // ---- Pipe system ----
-            int intakeCells  = 7;    // 100 mm stub  → ~14 mm/cell
+            int intakeCells    = 8;
-            int runnerCells  = 13;   // 180 mm runner → ~14 mm/cell
+            int reedCells      = 4;
-            int exhaustCells = 60;   // 850 mm total  → ~14 mm/cell
+            int transferCells  = 8;
+            int exhaustCells   = 60;
+            int totalCells = intakeCells + reedCells + transferCells + exhaustCells;
 
-            int totalCells = intakeCells + runnerCells + exhaustCells;
+            int[] pipeStart = {
-            int[] pipeStart = { 0,           intakeCells,              intakeCells + runnerCells };
+                0,
-            int[] pipeEnd   = { intakeCells, intakeCells + runnerCells, totalCells };
+                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];
 
-            // --- Intake ---
+            float intakeDia = 0.042f, reedDia = 0.040f, transferDia = 0.040f;
-            float intakeDia      = 0.042f;   // matches throttle body
+            intakePipeArea   = MathF.PI * 0.25f * intakeDia * intakeDia;
-            float intakeStubLen  = 0.100f;
+            reedPipeArea     = MathF.PI * 0.25f * reedDia * reedDia;
-            float intakeRunnerLen= 0.160f;   // shorter runner → less pumping loss
+            transferPipeArea = MathF.PI * 0.25f * transferDia * transferDia;
-            intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia;
 
             for (int i = 0; i < intakeCells; i++)
-            { area[i] = intakePipeArea; dx[i] = intakeStubLen / intakeCells; }
+            { area[i] = intakePipeArea; dx[i] = 0.100f / intakeCells; }
 
-            for (int i = intakeCells; i < intakeCells + runnerCells; i++)
+            for (int i = intakeCells; i < intakeCells + reedCells; i++)
-            { area[i] = intakePipeArea; dx[i] = intakeRunnerLen / runnerCells; }
+            { area[i] = reedPipeArea; dx[i] = 0.030f / reedCells; }
 
-            // Expansion chamber tuned for ~8 500 RPM power peak.
+            for (int i = intakeCells + reedCells; i < intakeCells + reedCells + transferCells; i++)
-            // Return-pulse travel distance = 0.5 × c_avg × (60 / RPM_target)
+            { area[i] = transferPipeArea; dx[i] = 0.200f / transferCells; }
-            // c_avg ≈ 480 m/s → distance = 0.5 × 480 × (60/8500) ≈ 1.69 m round-trip
-            // → one-way pipe length ≈ 0.84 m  (matches total below)
-            float headerDia    = 0.040f;  float headerLen    = 0.130f;   // shorter header → earlier pulse
-            float diffEndDia   = 0.070f;  float diffuserLen  = 0.250f;   // slightly narrower belly
-            float bellyDia     = 0.070f;  float bellyLen     = 0.220f;
-            float convEndDia   = 0.028f;  float convergentLen= 0.160f;   // longer convergent → stronger return pulse
-            float stingerDia   = 0.028f;  float stingerLen   = 0.080f;
-            // total = 0.13+0.25+0.22+0.16+0.08 = 0.84 m
 
-            exhaustHeaderArea  = MathF.PI * 0.25f * headerDia  * headerDia;
+            float hdrD = 0.040f, hdrL = 0.130f;
-            float bellyArea    = MathF.PI * 0.25f * bellyDia   * bellyDia;
+            float difEndD = 0.070f, difL = 0.250f;
-            float stingerArea  = MathF.PI * 0.25f * stingerDia * stingerDia;
+            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;
 
-            // Distribute cells proportionally by section length
+            int exhStart = intakeCells + reedCells + transferCells;
-            int headerCells     = Math.Max(1, (int)MathF.Round(exhaustCells * headerLen    / 0.84f));
+            int hdrC = (int)(exhaustCells * hdrL / totL);
-            int diffuserCells   = Math.Max(1, (int)MathF.Round(exhaustCells * diffuserLen  / 0.84f));
+            int difC = (int)(exhaustCells * difL / totL);
-            int bellyCells      = Math.Max(1, (int)MathF.Round(exhaustCells * bellyLen     / 0.84f));
+            int belC = (int)(exhaustCells * belL / totL);
-            int convergentCells = Math.Max(1, (int)MathF.Round(exhaustCells * convergentLen/ 0.84f));
+            int conC = (int)(exhaustCells * convL / totL);
-            int stingerCells    = exhaustCells - headerCells - diffuserCells
+            int stiC = exhaustCells - hdrC - difC - belC - conC;
-                                               - bellyCells  - convergentCells;
-            if (stingerCells < 1) stingerCells = 1;
 
-            int exhBase = intakeCells + runnerCells;
             int idx = 0;
-            for (int i = exhBase; i < totalCells; i++, idx++)
+            for (int i = exhStart; i < totalCells; i++)
             {
-                if (idx < headerCells)
+                if (idx < hdrC) { area[i] = exhaustHeaderArea; dx[i] = hdrL / hdrC; }
+                else if (idx < hdrC + difC)
                 {
-                    area[i] = exhaustHeaderArea;
+                    float t = (idx - hdrC) / (float)(difC - 1);
-                    dx[i]   = headerLen / headerCells;
+                    float dia = hdrD + (difEndD - hdrD) * t;
+                    area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = difL / difC;
                 }
-                else if (idx < headerCells + diffuserCells)
+                else if (idx < hdrC + difC + belC) { area[i] = bellyArea; dx[i] = belL / belC; }
+                else if (idx < hdrC + difC + belC + conC)
                 {
-                    float t  = (idx - headerCells) / (float)(diffuserCells - 1);
+                    float t = (idx - hdrC - difC - belC) / (float)(conC - 1);
-                    // Smooth cosine taper instead of linear for better wave reflection
+                    float dia = difEndD + (convEndD - difEndD) * t;
-                    float ct = 0.5f * (1f - MathF.Cos(MathF.PI * t));
+                    area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = convL / conC;
-                    float dia = headerDia + (diffEndDia - headerDia) * ct;
-                    area[i] = MathF.PI * 0.25f * dia * dia;
-                    dx[i]   = diffuserLen / diffuserCells;
-                }
-                else if (idx < headerCells + diffuserCells + bellyCells)
-                {
-                    area[i] = bellyArea;
-                    dx[i]   = bellyLen / bellyCells;
-                }
-                else if (idx < headerCells + diffuserCells + bellyCells + convergentCells)
-                {
-                    float t  = (idx - headerCells - diffuserCells - bellyCells)
-                               / (float)(convergentCells - 1);
-                    // Steeper cosine convergent for a sharper return pulse
-                    float ct = 0.5f * (1f - MathF.Cos(MathF.PI * t));
-                    float dia = bellyDia + (convEndDia - bellyDia) * ct;
-                    area[i] = MathF.PI * 0.25f * dia * dia;
-                    dx[i]   = convergentLen / convergentCells;
-                }
-                else
-                {
-                    area[i] = stingerArea;
-                    dx[i]   = stingerLen / stingerCells;
                 }
+                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;   // slightly less damping → stronger pulses
+            pipeSystem.DampingMultiplier = 0.8f;
             pipeSystem.EnergyRelaxationRate = 0.4f;
-            pipeSystem.AmbientPressure      = 101325f;
 
-            // ── 0-D Volumes ──────────────────────────────────────────────────────
+            // ---- Volumes ----
-            // Intake plenum: acts as a small airbox resonator (8 cc)
+            intakePlenum = new Volume0D(0.5e-3f, 101325f, 300f);
-            intakePlenum = new Volume0D(8e-3f, 101325f, 300f);
             plenumInlet  = intakePlenum.CreatePort();
             plenumOutlet = intakePlenum.CreatePort();
 
-            // Exhaust silencer volume: 600 cc is realistic for a small-bore muffler
+            // ---- Boundary system ----
-            exhaustMuffler = new Volume0D(600e-6f, 101325f, 650f);
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 6, maxOpenEnds: 2);
-            mufflerIn  = exhaustMuffler.CreatePort();
+            throttleAreaIdx   = 0;
-            mufflerOut = exhaustMuffler.CreatePort();
+            reedInletIdx      = 1;
+            reedOutletIdx     = 2;
+            transferInletIdx  = 3;
+            transferOutletIdx = 4;
+            exhaustValveIdx   = 5;
 
-            // ── Boundary system ───────────────────────────────────────────────────
+            boundaries.AddOpenEnd(0, true,  101325f, intakePipeArea);
-            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
+            intakeOpenIdx = 0;
-            throttleAreaIdx  = 0;
+            boundaries.AddOpenEnd(3, false, 101325f, stingerArea);
-            plenumRunnerIdx  = 1;
-            intakeValveIdx   = 2;
-            exhaustValveIdx  = 3;
-
-            // Open ends: atmosphere at both extremes
-            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true,  101325f, intakePipeArea);
-            intakeOpenIdx  = 0;
-            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, stingerArea);
             exhaustOpenIdx = 1;
 
-            // Orifices: throttle → plenum → runner → cylinder → exhaust pipe
+            boundaries.AddOrifice(plenumInlet,           0, false, throttleAreaIdx, 0.72f);
-            boundaries.AddOrifice(plenumInlet,          0, false, throttleAreaIdx, 0.72f);
+            boundaries.AddOrifice(plenumOutlet,          1, true,  reedInletIdx,   1.0f);
-            boundaries.AddOrifice(plenumOutlet,         1, true,  plenumRunnerIdx, 1.00f);
+            boundaries.AddOrifice(crankcase.IntakePort,  1, false, reedOutletIdx,  0.9f);
-            boundaries.AddOrifice(cylinder.IntakePort,  1, false, intakeValveIdx,  0.68f);
+            boundaries.AddOrifice(crankcase.TransferPort,2, true,  transferInletIdx,1.0f);
-            boundaries.AddOrifice(cylinder.ExhaustPort, 2, true,  exhaustValveIdx, 0.70f);
+            boundaries.AddOrifice(cylinder.IntakePort,   2, false, transferOutletIdx,1.0f);
+            boundaries.AddOrifice(cylinder.ExhaustPort,  3, true,  exhaustValveIdx, 0.7f);
 
-            orificeAreas = new float[4];
+            orificeAreas = new float[6];
-            orificeAreas[plenumRunnerIdx] = intakePipeArea;  // runner always fully open
+            orificeAreas[reedInletIdx]      = reedPipeArea;
+            orificeAreas[reedOutletIdx]     = 0f;
+            orificeAreas[transferInletIdx]  = transferPipeArea;
+            orificeAreas[transferOutletIdx] = 0f;
 
-            // ── Solver ────────────────────────────────────────────────────────────
+            // ---- Solver ----
-            // SubStepCount = 4 keeps CFL ≤ 1 for 5 mm cells at 44 100 Hz
+            solver = new Solver { SubStepCount = 4 };
-            solver = new Solver { SubStepCount = 4, EnableProfiling = false };
             solver.SetTimeStep(dt);
             solver.SetPipeSystem(pipeSystem);
             solver.SetBoundarySystem(boundaries);
             solver.AddComponent(cylinder);
+            solver.AddComponent(crankcase);
             solver.AddComponent(intakePlenum);
-            solver.AddComponent(exhaustMuffler);
 
-            // ── Sound ─────────────────────────────────────────────────────────────
+            // ---- 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("125cc Two-Stroke – expansion chamber tuned for ~8 500 RPM power peak");
+            Console.WriteLine("Two‑Stroke engine ready.");
-            Console.WriteLine($"  Exhaust cells: {exhaustCells}  |  header {headerCells}  diffuser {diffuserCells}" +
-                              $"  belly {bellyCells}  convergent {convergentCells}  stinger {stingerCells}");
         }
 
         public override float Process()
         {
-            float engineRpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            const float reedMargin = 200f;
+            if (crankcase.Pressure < intakePlenum.Pressure - reedMargin)
+                reedOpen = true;
+            else if (crankcase.Pressure > intakePlenum.Pressure + reedMargin)
+                reedOpen = false;
 
-            vehicle.ClutchInput = Clutch;
+            float throttledFraction = Throttle;
+            if (throttledFraction < 0.001f) throttledFraction = 0f;
+            throttledFraction = Math.Clamp(throttledFraction, 0f, 1f);
+            float throttledArea = maxThrottleArea * throttledFraction;
 
-            var (clutchTorque, effectiveInertia) = vehicle.Update(engineRpm, crankshaft.Inertia, (float)dt);
+            orificeAreas[throttleAreaIdx] = throttledArea;
-            crankshaft.SetEffectiveInertia(effectiveInertia);
+            orificeAreas[reedOutletIdx] = reedOpen ? reedPipeArea : 0f;
-            crankshaft.SetLoadTorque(clutchTorque);
+            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);
-            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
-            orificeAreas[throttleAreaIdx]  = throttledArea;
-            orificeAreas[intakeValveIdx]   = cylinder.IntakeValveArea;
-            orificeAreas[exhaustValveIdx]  = cylinder.ExhaustValveArea;
-            boundaries.SetOrificeAreas(orificeAreas);
 
             solver.Step();
             stepCount++;
@@ -273,76 +237,50 @@ namespace FluidSim.Tests
 
             if (stepCount % 2000 == 0)
             {
-                float rpm     = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+                float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-                float powerKw = crankshaft.AveragePower * 1e-3f;
+                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}");
-                float torqueNm = crankshaft.AverageTorque;
-                Console.WriteLine($"Step {stepCount,7} | RPM={rpm,6:F0} | Power={powerKw,5:F2} kW" +
-                                  $" | Torque={torqueNm,5:F1} Nm | Gear={vehicle.CurrentGear}" +
-                                  $" | Speed={vehicle.SpeedKmh,4:F0} km/h");
             }
 
             return reverb.Process((intakeDry + exhaustDry) * 0.5f);
         }
 
-        // ── Drawing ───────────────────────────────────────────────────────────────
         public override void Draw(RenderWindow target)
         {
             float winW = target.GetView().Size.X;
             float winH = target.GetView().Size.Y;
 
-            float intakeY  = winH / 2f - 40f;
+            float startX = 40f;
-            float exhaustY = winH / 2f + 80f;
+            float endX = winW - 80f;
-            float openEndX = 40f;
 
-            // Intake stub
+            DrawPipe(target, pipeSystem, 0, winH * 0.25f, startX, startX + 120f);
-            float x = openEndX;
-            float w = 120f;
-            DrawPipe(target, pipeSystem, 0, intakeY, x, x + w);
-
-            // Throttle body
-            float throttleX = x + w + 5f;
             var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
             {
                 FillColor = Color.Yellow,
-                Position  = new Vector2f(throttleX, intakeY - 15f)
+                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);
 
-            // Plenum
+            float reedStartX = plenX + 70f;
-            float plenW = 40f, plenH = 60f;
+            DrawPipe(target, pipeSystem, 1, winH * 0.25f, reedStartX, reedStartX + 30f);
-            float plenX = throttleX + 10f;
-            DrawVolume(target, intakePlenum, plenX + plenW / 2f, intakeY - plenH / 2f, plenW, plenH);
 
-            // Runner
+            float transStartX = reedStartX + 40f;
-            float runnerStartX = plenX + plenW + 5f;
+            DrawPipe(target, pipeSystem, 2, winH * 0.45f, transStartX, transStartX + 120f);
-            DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerStartX + 100f);
 
-            // Cylinder
+            float cylCX = transStartX + 180f;
-            float cylCX  = runnerStartX + 150f;
+            float cylTopY = winH * 0.45f - 90f;
-            float cylTopY = intakeY - 120f;
             DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f);
 
-            // Exhaust pipe (expansion chamber)
+            float exhStartX = cylCX + 60f;
-            float exhStartX = cylCX + 40f + 20f;
+            DrawPipe(target, pipeSystem, 3, winH * 0.65f, exhStartX, endX, areaScale: 800f);
-            DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, winW - 60f, areaScale: 800f);
 
-            // HUD labels
+            float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-            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;
-
-            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);
-            DrawLabel(target, $"Torque: {torqueNm:F1} Nm",new Vector2f(20, 140), Color.White, 20);
-
-            string gearText = vehicle.CurrentGear == 0 ? "N" : vehicle.CurrentGear.ToString();
-            DrawLabel(target, $"Gear: {gearText}",        new Vector2f(20, 162), Color.Cyan, 20);
-            DrawLabel(target, $"Speed: {vehicle.SpeedKmh:F0} km/h",
-                                                          new Vector2f(20, 184), Color.Cyan, 20);
-            DrawLabel(target, vehicle.Engagement > 0.99f ? "Clutch: Locked" : "Clutch: Slipping",
-                                                          new Vector2f(20, 204), Color.Cyan, 14);
-
-            // Dyno curve
             UpdateDynoCurve(rpm, powerKw, torqueNm);
             DrawDynoCurve(target, winW - 410f, winH - 260f, 400f, 250f, rpm, powerKw);
         }