From 16498f8041f198563d2844011974d6fee9a88fdf Mon Sep 17 00:00:00 2001 From: max Date: Wed, 10 Jun 2026 00:13:14 +0200 Subject: [PATCH] crank case based two stroke testing --- Components/Crankcase.cs | 150 +++++++++++++ Components/Crankshaft.cs | 1 + Components/TwoStrokeCylinder.cs | 199 ++++++++--------- Scenarios/TwoStrokeScenario.cs | 364 +++++++++++++------------------- 4 files changed, 406 insertions(+), 308 deletions(-) create mode 100644 Components/Crankcase.cs diff --git a/Components/Crankcase.cs b/Components/Crankcase.cs new file mode 100644 index 0000000..05edb1b --- /dev/null +++ b/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 _ports; + public IReadOnlyList 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 { 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; + } + } + } + } +} \ No newline at end of file diff --git a/Components/Crankshaft.cs b/Components/Crankshaft.cs index e136c92..8c393b7 100644 --- a/Components/Crankshaft.cs +++ b/Components/Crankshaft.cs @@ -27,6 +27,7 @@ namespace FluidSim.Components /// Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke. public float CycleLength { get; set; } = 4f * MathF.PI; + public float CrankAngleRad => CrankAngle; public Crankshaft(float initialRPM = 400f) { diff --git a/Components/TwoStrokeCylinder.cs b/Components/TwoStrokeCylinder.cs index 5c928bc..3a4bf8c 100644 --- a/Components/TwoStrokeCylinder.cs +++ b/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 { /// - /// Two-stroke cylinder with symmetrical port timings centred on BDC (180°). - /// - /// 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. + /// Two‑stroke cylinder with forced symmetrical port timings around BDC (180°). + /// Uses crankcase back‑pressure for accurate pumping work. /// public class TwoStrokeCylinder : EngineCylinder { - // ── Port timing read-outs (degrees, 0 = TDC) ─────────────────────────── - public float IVO => 180f - TransferDuration / 2f; // transfer opens - public float IVC => 180f + TransferDuration / 2f; // transfer closes - public float EVO => 180f - ExhaustDuration / 2f; // exhaust opens - public float EVC => 180f + ExhaustDuration / 2f; // exhaust closes + // --- 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; - // ── Configurable durations ────────────────────────────────────────────── - public float TransferDuration { get; } // default: 155° - public float ExhaustDuration { get; } // default: 195° + private readonly float transferDuration; // degrees + private readonly float exhaustDuration; // degrees - // Fraction of port-open duration used for ramp-up / ramp-down. - // 0.15 → port at full area for the middle 70 % of open time. - private const float RampFraction = 0.15f; + // --- Crankcase reference --- + private Crankcase? _crankcase; protected override float CycleLengthRad => 2f * MathF.PI; - protected override float MaxCycleDeg => 360f; + protected override float MaxCycleDeg => 360f; public override float IntakeValveArea => - MathF.PI * IntakeValveDiameter - * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift); - + MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift); public override float ExhaustValveArea => - MathF.PI * ExhaustValveDiameter - * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift); + MathF.PI * ExhaustValveDiameter * 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; - ExhaustDuration = exhaustDuration; + this.transferDuration = transferDuration; + this.exhaustDuration = exhaustDuration; if (EVO >= IVO) - throw new ArgumentException( - $"Exhaust must open before transfer port. " + - $"EVO={EVO:F1}° must be less than IVO={IVO:F1}°. " + - $"Increase exhaustDuration or decrease transferDuration."); + throw new ArgumentException("Exhaust must open before transfer port."); } - // ── Valve lift profile ────────────────────────────────────────────────── - /// - /// Smooth trapezoidal lift: fast ramp (15 % of duration), flat top (70 %), - /// fast ramp-down (15 %). Ramps use a smoothstep (3t²-2t³) curve so the - /// area derivative is C1-continuous (no kink at ramp/plateau boundaries). - /// - private static float ValveLift(float thetaDeg, float opens, float closes, float peakLift) + public void SetCrankcase(Crankcase crankcase) + { + _crankcase = crankcase; + } + + // ----- Valve lift ----- + 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; - // a port that crosses 360° would need closes+360). - float effectiveClose = closes < opens ? closes + 360f : closes; - float duration = effectiveClose - opens; + float effectiveOpen = opens; + float effectiveClose = closes; + if (closes < opens) effectiveClose += 360f; + float duration = effectiveClose - effectiveOpen; if (duration <= 0f) return 0f; - // Map deg into the same number-line as opens/effectiveClose - float mapped = deg < opens ? deg + 360f : deg; + float mapped = deg; + if (mapped < opens) mapped += 360f; if (mapped < opens || mapped > effectiveClose) return 0f; - float rampDur = duration * RampFraction; - float holdEnd = effectiveClose - rampDur; + float rampDur = duration * 0.25f; + 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 - holdEnd) / rampDur; + float t = (mapped - (opens + rampDur + holdDur)) / 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) ────────────────── - // At IVC the cylinder is sealed; whatever air is trapped is what we burn. - if (CrossedAngle(prevDeg, currDeg, IVC)) + // Transfer port closing → fuel injection + if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC) { trappedAirMass = _airMass; - fuelMass = trappedAirMass / StoichiometricAFR; - fuelInjected = true; + fuelMass = trappedAirMass / StoichiometricAFR; + fuelInjected = true; } - // ── Ignition ─────────────────────────────────────────────────────── - // SparkAdvance default is ~22° BTDC on the base class; scenario can override. - float sparkAngle = (360f - SparkAdvance) % 360f; + // 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 (CrossedAngle(prevDeg, currDeg, sparkAngle) && !combustionActive && fuelInjected) + if (crossedSpark && !combustionActive && fuelInjected) { if (_random.NextDouble() < MisfireProbability) { @@ -126,58 +113,80 @@ namespace FluidSim.Components } else { - combustionActive = true; - burnFraction = 0f; + combustionActive = true; 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 - || angleSinceSpark > WiebeDuration + WiebeStart + SparkAdvance; - - if (burnComplete) + if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance)) { - newFraction = 1f; - combustionActive = false; - fuelInjected = false; - float totalMass = _airMass + _exhaustMass; - _airMass = 0f; - _exhaustMass = totalMass; + 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; + _exhaustMass += fuelMass * dFraction; + burnFraction = newFraction; } } } - // ── Helper: did the crank cross a target angle this step? ─────────────── - /// - /// Returns true if the crank swept through going - /// from to in a single step. - /// Handles wrap-around at 360°. - /// - private static bool CrossedAngle(float prev, float curr, float target) + // ----- Override torque calculation to use crankcase back‑pressure ----- + public new void PreStep(float dt) { - // Normal case (no wrap) - if (curr >= prev) - return prev < target && target <= curr; + // Speed‑dependent spark advance + float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI); + SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f); - // Wrapped past 360° → two intervals to check - return prev < target || target <= curr; + 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; } } } \ No newline at end of file diff --git a/Scenarios/TwoStrokeScenario.cs b/Scenarios/TwoStrokeScenario.cs index f450946..b63fb1c 100644 --- a/Scenarios/TwoStrokeScenario.cs +++ b/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, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx; + private int throttleAreaIdx, reedInletIdx, reedOutletIdx, + transferInletIdx, transferOutletIdx, 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 intakePipeArea, exhaustHeaderArea; - - public override void ShiftUp() => vehicle.ShiftUp(); - public override void ShiftDown() => vehicle.ShiftDown(); + private float maxThrottleArea; + private float intakePipeArea, reedPipeArea, transferPipeArea, exhaustHeaderArea; + private bool reedOpen; public override void Initialize(int sampleRate) { dt = 1.0 / sampleRate; - // ── Vehicle ────────────────────────────────────────────────────────── - vehicle = new Vehicle(); + maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm); - // ── Throttle body: 42 mm – wider to reduce high-RPM intake restriction ── - _maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm); - - // ── Crankshaft ─────────────────────────────────────────────────────── - // Lighter flywheel for quicker revving; friction tuned to ~0.5 kW loss at idle - crankshaft = new Crankshaft(2000); - 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; + // ---- 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, - transferDuration, exhaustDuration, - crankshaft) + transferDur, exhaustDur, crankshaft) { - IntakeValveDiameter = 0.042f, // matched to intake pipe - IntakeValveLift = 0.015f, + // FIX: realistic transfer port diameter (was 40mm) + IntakeValveDiameter = 0.030f, // 30 mm + IntakeValveLift = 0.010f, ExhaustValveDiameter = 0.040f, - ExhaustValveLift = 0.013f + ExhaustValveLift = 0.010f }; - // ── Pipe geometry ──────────────────────────────────────────────────── - // - // Layout (all lengths in mm): - // Intake path: airbox stub 100 mm | runner 180 mm - // Exhaust path: expansion chamber tuned to ~9 000 RPM power peak - // header 170 mm Ø 40 mm - // 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 + // ---- Crankcase ---- + float crankRadius = stroke * 0.5f; + float ccClearance = 150e-6f; + crankcase = new Crankcase(crankshaft, crankRadius, conRod, bore, + ccClearance, 101325f, 300f); + cylinder.SetCrankcase(crankcase); - // --- Cell counts --- - int intakeCells = 7; // 100 mm stub → ~14 mm/cell - int runnerCells = 13; // 180 mm runner → ~14 mm/cell - int exhaustCells = 60; // 850 mm total → ~14 mm/cell + // ---- Pipe system ---- + int intakeCells = 8; + int reedCells = 4; + int transferCells = 8; + int exhaustCells = 60; + int totalCells = intakeCells + reedCells + transferCells + exhaustCells; - int totalCells = intakeCells + runnerCells + exhaustCells; - int[] pipeStart = { 0, intakeCells, intakeCells + runnerCells }; - int[] pipeEnd = { intakeCells, intakeCells + runnerCells, totalCells }; + 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]; - // --- Intake --- - float intakeDia = 0.042f; // matches throttle body - float intakeStubLen = 0.100f; - float intakeRunnerLen= 0.160f; // shorter runner → less pumping loss - intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia; + 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] = intakeStubLen / intakeCells; } + { area[i] = intakePipeArea; dx[i] = 0.100f / intakeCells; } - for (int i = intakeCells; i < intakeCells + runnerCells; i++) - { area[i] = intakePipeArea; dx[i] = intakeRunnerLen / runnerCells; } + for (int i = intakeCells; i < intakeCells + reedCells; i++) + { area[i] = reedPipeArea; dx[i] = 0.030f / reedCells; } - // Expansion chamber tuned for ~8 500 RPM power peak. - // Return-pulse travel distance = 0.5 × c_avg × (60 / RPM_target) - // 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 + for (int i = intakeCells + reedCells; i < intakeCells + reedCells + transferCells; i++) + { area[i] = transferPipeArea; dx[i] = 0.200f / transferCells; } - exhaustHeaderArea = MathF.PI * 0.25f * headerDia * headerDia; - float bellyArea = MathF.PI * 0.25f * bellyDia * bellyDia; - float stingerArea = MathF.PI * 0.25f * stingerDia * stingerDia; + 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; - // Distribute cells proportionally by section length - int headerCells = Math.Max(1, (int)MathF.Round(exhaustCells * headerLen / 0.84f)); - int diffuserCells = Math.Max(1, (int)MathF.Round(exhaustCells * diffuserLen / 0.84f)); - int bellyCells = Math.Max(1, (int)MathF.Round(exhaustCells * bellyLen / 0.84f)); - int convergentCells = Math.Max(1, (int)MathF.Round(exhaustCells * convergentLen/ 0.84f)); - int stingerCells = exhaustCells - headerCells - diffuserCells - - bellyCells - convergentCells; - if (stingerCells < 1) stingerCells = 1; + 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 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; - dx[i] = headerLen / headerCells; + 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 < 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); - // Smooth cosine taper instead of linear for better wave reflection - float ct = 0.5f * (1f - MathF.Cos(MathF.PI * t)); - 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; + 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; // slightly less damping → stronger pulses + pipeSystem.DampingMultiplier = 0.8f; pipeSystem.EnergyRelaxationRate = 0.4f; - pipeSystem.AmbientPressure = 101325f; - // ── 0-D Volumes ────────────────────────────────────────────────────── - // Intake plenum: acts as a small airbox resonator (8 cc) - intakePlenum = new Volume0D(8e-3f, 101325f, 300f); + // ---- Volumes ---- + intakePlenum = new Volume0D(0.5e-3f, 101325f, 300f); plenumInlet = intakePlenum.CreatePort(); plenumOutlet = intakePlenum.CreatePort(); - // Exhaust silencer volume: 600 cc is realistic for a small-bore muffler - exhaustMuffler = new Volume0D(600e-6f, 101325f, 650f); - mufflerIn = exhaustMuffler.CreatePort(); - mufflerOut = exhaustMuffler.CreatePort(); + // ---- Boundary system ---- + boundaries = new BoundarySystem(pipeSystem, maxOrifices: 6, maxOpenEnds: 2); + throttleAreaIdx = 0; + reedInletIdx = 1; + reedOutletIdx = 2; + transferInletIdx = 3; + transferOutletIdx = 4; + exhaustValveIdx = 5; - // ── Boundary system ─────────────────────────────────────────────────── - boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2); - throttleAreaIdx = 0; - 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); + boundaries.AddOpenEnd(0, true, 101325f, intakePipeArea); + intakeOpenIdx = 0; + boundaries.AddOpenEnd(3, false, 101325f, stingerArea); exhaustOpenIdx = 1; - // Orifices: throttle → plenum → runner → cylinder → exhaust pipe - boundaries.AddOrifice(plenumInlet, 0, false, throttleAreaIdx, 0.72f); - boundaries.AddOrifice(plenumOutlet, 1, true, plenumRunnerIdx, 1.00f); - boundaries.AddOrifice(cylinder.IntakePort, 1, false, intakeValveIdx, 0.68f); - boundaries.AddOrifice(cylinder.ExhaustPort, 2, true, exhaustValveIdx, 0.70f); + 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[4]; - orificeAreas[plenumRunnerIdx] = intakePipeArea; // runner always fully open + orificeAreas = new float[6]; + orificeAreas[reedInletIdx] = reedPipeArea; + orificeAreas[reedOutletIdx] = 0f; + orificeAreas[transferInletIdx] = transferPipeArea; + orificeAreas[transferOutletIdx] = 0f; - // ── Solver ──────────────────────────────────────────────────────────── - // SubStepCount = 4 keeps CFL ≤ 1 for 5 mm cells at 44 100 Hz - solver = new Solver { SubStepCount = 4, EnableProfiling = false }; + // ---- Solver ---- + solver = new Solver { SubStepCount = 4 }; 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($" Exhaust cells: {exhaustCells} | header {headerCells} diffuser {diffuserCells}" + - $" belly {bellyCells} convergent {convergentCells} stinger {stingerCells}"); + Console.WriteLine("Two‑Stroke engine ready."); } 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); - crankshaft.SetEffectiveInertia(effectiveInertia); - crankshaft.SetLoadTorque(clutchTorque); + 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); - - float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f); - orificeAreas[throttleAreaIdx] = throttledArea; - orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea; - orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea; - boundaries.SetOrificeAreas(orificeAreas); + crankcase.PreStep((float)dt); solver.Step(); stepCount++; @@ -273,76 +237,50 @@ namespace FluidSim.Tests if (stepCount % 2000 == 0) { - float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI); - float powerKw = crankshaft.AveragePower * 1e-3f; - 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"); + 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); } - // ── 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 exhaustY = winH / 2f + 80f; - float openEndX = 40f; + float startX = 40f; + float endX = winW - 80f; - // Intake stub - float x = openEndX; - float w = 120f; - DrawPipe(target, pipeSystem, 0, intakeY, x, x + w); - - // Throttle body - float throttleX = x + w + 5f; + DrawPipe(target, pipeSystem, 0, winH * 0.25f, startX, startX + 120f); 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 plenW = 40f, plenH = 60f; - float plenX = throttleX + 10f; - DrawVolume(target, intakePlenum, plenX + plenW / 2f, intakeY - plenH / 2f, plenW, plenH); + float reedStartX = plenX + 70f; + DrawPipe(target, pipeSystem, 1, winH * 0.25f, reedStartX, reedStartX + 30f); - // Runner - float runnerStartX = plenX + plenW + 5f; - DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerStartX + 100f); + float transStartX = reedStartX + 40f; + DrawPipe(target, pipeSystem, 2, winH * 0.45f, transStartX, transStartX + 120f); - // Cylinder - float cylCX = runnerStartX + 150f; - float cylTopY = intakeY - 120f; + float cylCX = transStartX + 180f; + float cylTopY = winH * 0.45f - 90f; DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f); - // Exhaust pipe (expansion chamber) - float exhStartX = cylCX + 40f + 20f; - DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, winW - 60f, areaScale: 800f); + float exhStartX = cylCX + 60f; + DrawPipe(target, pipeSystem, 3, winH * 0.65f, exhStartX, endX, 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); }