diff --git a/Components/Crankshaft.cs b/Components/Crankshaft.cs
index 0af1361..e136c92 100644
--- a/Components/Crankshaft.cs
+++ b/Components/Crankshaft.cs
@@ -4,33 +4,30 @@ namespace FluidSim.Components
{
public class Crankshaft
{
- public float AngularVelocity; // rad/s
- public float CrankAngle; // rad, 0 … 4π
+ public float AngularVelocity;
+ public float CrankAngle;
public float PreviousAngle;
- public float Inertia = 0.2f; // kg·m²
- public float FrictionConstant; // N·m
- public float FrictionViscous; // N·m per rad/s
-
+ public float Inertia = 0.2f;
+ public float FrictionConstant;
+ public float FrictionViscous;
public float LastNetTorque { get; private set; }
- public float AveragePower { get; private set; } // smoothed, watts
- public float AverageTorque { get; private set; } // smoothed, Nm
+ public float AveragePower { get; private set; }
+ 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 _powerBufCount;
+ private int _powerBufIdx, _powerBufCount;
private float _powerBufSum;
-
- // Torque averaging buffer (same size as power buffer)
private readonly float[] _torqueBuffer;
- private int _torqueBufIdx;
- private int _torqueBufCount;
+ private int _torqueBufIdx, _torqueBufCount;
private float _torqueBufSum;
+ /// Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke.
+ public float CycleLength { get; set; } = 4f * MathF.PI;
+
public Crankshaft(float initialRPM = 400f)
{
AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
@@ -43,9 +40,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 +58,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)
- CrankAngle -= 4f * MathF.PI;
+ if (CrankAngle >= CycleLength)
+ 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;
diff --git a/Components/Cylinder.cs b/Components/Cylinder.cs
index 8b18342..6dc57cd 100644
--- a/Components/Cylinder.cs
+++ b/Components/Cylinder.cs
@@ -1,99 +1,25 @@
using System;
-using System.Collections.Generic;
-using FluidSim.Interfaces;
+using FluidSim.Components; // if needed
namespace FluidSim.Components
{
- public class Cylinder : IComponent
+ public class Cylinder : EngineCylinder
{
- public Port IntakePort { get; }
- public Port ExhaustPort { get; }
- public Crankshaft Crankshaft { get; }
+ public float IVO, IVC, EVO, EVC; // degrees in a 720° cycle
- private readonly Port[] _ports;
- IReadOnlyList IComponent.Ports => _ports;
+ protected override float CycleLengthRad => 4f * MathF.PI;
+ protected override float MaxCycleDeg => 720f;
- public float Bore { get; }
- public float Stroke { get; }
- public float ConRodLength { get; }
- public float CompressionRatio { get; }
-
- 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 override float IntakeValveArea =>
+ MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
+ public override float ExhaustValveArea =>
+ MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
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;
- }
- duration = effectiveClose - effectiveOpen;
+ if (closes < opens) effectiveClose += 720f;
+ float duration = effectiveClose - effectiveOpen;
if (duration <= 0f) return 0f;
float mapped = deg;
@@ -136,43 +57,9 @@ namespace FluidSim.Components
return 0f;
}
- public float IntakeValveArea =>
- MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
- public float ExhaustValveArea =>
- MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-
- private float Wiebe(float angleSinceSpark)
+ protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
{
- if (angleSinceSpark < WiebeStart) return 0f;
- 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
+ // Intake closing → fuel injection
if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
{
trappedAirMass = _airMass;
@@ -180,11 +67,14 @@ namespace FluidSim.Components
fuelInjected = true;
}
- // Spark
- float sparkAngle = 0f - SparkAdvance;
- if (sparkAngle < 0f) sparkAngle += 720f;
- bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) ||
- (prevDeg > sparkAngle + 360f && currDeg < sparkAngle);
+ // Spark – occurs at TDC (0°) minus advance, every 720°
+ float sparkAngle = (0f - SparkAdvance + 720f) % 720f;
+ 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)
@@ -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;
}
}
}
\ No newline at end of file
diff --git a/Components/EngineCylinder.cs b/Components/EngineCylinder.cs
new file mode 100644
index 0000000..508e88c
--- /dev/null
+++ b/Components/EngineCylinder.cs
@@ -0,0 +1,203 @@
+using System;
+using System.Collections.Generic;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+ /// Common base for all reciprocating engine cylinders.
+ public abstract class EngineCylinder : IComponent
+ {
+ public Port IntakePort { get; }
+ public Port ExhaustPort { get; }
+ public Crankshaft Crankshaft { get; }
+
+ private readonly Port[] _ports;
+ IReadOnlyList 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;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/TwoStrokeCylinder.cs b/Components/TwoStrokeCylinder.cs
new file mode 100644
index 0000000..ae6d9d3
--- /dev/null
+++ b/Components/TwoStrokeCylinder.cs
@@ -0,0 +1,140 @@
+using System;
+
+namespace FluidSim.Components
+{
+ ///
+ /// Two‑stroke cylinder with forced symmetrical port timings around BDC (180°).
+ /// All angles are in degrees within a 360° cycle.
+ ///
+ public class TwoStrokeCylinder : EngineCylinder
+ {
+ // --- Public read‑only properties for drawing ---
+ 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 (set in constructor) ---
+ private readonly float transferDuration; // e.g. 120°
+ private readonly float exhaustDuration; // e.g. 180°
+
+ 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);
+
+ ///
+ /// Create a two‑stroke cylinder with forced symmetrical port timing.
+ ///
+ /// Total transfer port open duration in degrees (e.g. 120°).
+ /// Total exhaust port open duration in degrees (e.g. 180°).
+ 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;
+
+ // Safety check: exhaust must open before transfer
+ if (EVO >= IVO)
+ throw new ArgumentException("Exhaust must open before transfer port (exhaust duration > transfer duration).");
+ }
+
+ // ----- Valve lift – same implementation, now uses the computed IVO/IVC/EVO/EVC -----
+ 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;
+ }
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/Vehicle.cs b/Components/Vehicle.cs
new file mode 100644
index 0000000..f7c0962
--- /dev/null
+++ b/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;
+
+ /// Time constant for clutch engagement smoothing (seconds).
+ 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 ----
+ /// Peak clutch friction torque (Nm) when fully engaged at high RPM.
+ public float BaseMaxTorque = 80f; // much lower than before
+
+ /// Stiffness when slipping (Nm per rad/s). Lower = softer engagement.
+ public float ClutchStiffness = 50f; // very soft
+
+ /// Below this engine RPM, the clutch torque is progressively reduced to prevent stalling.
+ 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;
+ }
+}
\ No newline at end of file
diff --git a/Program.cs b/Program.cs
index c2668e2..a02a562 100644
--- a/Program.cs
+++ b/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;
}
}
\ No newline at end of file
diff --git a/Scenarios/Scenario.cs b/Scenarios/Scenario.cs
index 7c170de..b046103 100644
--- a/Scenarios/Scenario.cs
+++ b/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);
- float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
- radii[i] = baseRadius * (1f + dev * 2f);
- if (radii[i] < 2f) radii[i] = 2f;
+
+ 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;
}
diff --git a/Scenarios/TwoStrokeScenario.cs b/Scenarios/TwoStrokeScenario.cs
new file mode 100644
index 0000000..1ed0ca2
--- /dev/null
+++ b/Scenarios/TwoStrokeScenario.cs
@@ -0,0 +1,293 @@
+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 PipeSystem pipeSystem;
+ private BoundarySystem boundaries;
+ private Solver solver;
+
+ private Volume0D intakePlenum;
+ private Port plenumInlet, plenumOutlet;
+ private Volume0D exhaustMuffler;
+ private Port mufflerIn, mufflerOut;
+
+ private Vehicle vehicle;
+
+ private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, 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, exhaustHeaderArea;
+
+ // -- Override shift from Scenario base class --
+ public override void ShiftUp() => vehicle.ShiftUp();
+ public override void ShiftDown() => vehicle.ShiftDown();
+
+ public override void Initialize(int sampleRate)
+ {
+ dt = 1.0 / sampleRate;
+
+ // ---- Vehicle ----
+ vehicle = new Vehicle();
+
+ // ---- Throttle (38 mm) ----
+ _maxThrottleArea = (float)Units.AreaFromDiameter(38 * Units.mm);
+
+ // ---- Crankshaft ----
+ crankshaft = new Crankshaft(2000);
+ crankshaft.CycleLength = 2f * MathF.PI; // two‑stroke
+ crankshaft.Inertia = 0.05f; // engine's own inertia (light)
+ crankshaft.FrictionConstant = 2.5f;
+ crankshaft.FrictionViscous = 0.0015f;
+
+ // ---- Cylinder (125cc) ----
+ float bore = 0.054f, stroke = 0.0545f, conRod = 0.109f, compRatio = 12.5f;
+
+ // Symmetric durations (around BDC)
+ float transferDuration = 130f; // 130°
+ float exhaustDuration = 190f; // 190°
+
+ cylinder = new TwoStrokeCylinder(bore, stroke, conRod, compRatio,
+ transferDuration, exhaustDuration,
+ crankshaft)
+ {
+ IntakeValveDiameter = 0.038f,
+ IntakeValveLift = 0.010f,
+ ExhaustValveDiameter = 0.040f,
+ ExhaustValveLift = 0.010f
+ };
+
+ // ---- Pipe system (60 exhaust cells, simple diffuser) ----
+ int intakeCells = 8;
+ int runnerCells = 8;
+ int exhaustCells = 60;
+ int totalCells = intakeCells + runnerCells + exhaustCells;
+ int[] pipeStart = { 0, intakeCells, intakeCells + runnerCells };
+ int[] pipeEnd = { intakeCells, intakeCells + runnerCells, totalCells };
+
+ float[] area = new float[totalCells];
+ float[] dx = new float[totalCells];
+
+ float intakeDia = 0.038f;
+ float intakeLenBefore = 0.15f;
+ float intakeLenRunner = 0.20f;
+ intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia;
+
+ // Single‑stage diffuser – 840 mm total, easy to tune
+ float headerDia = 0.042f, headerLen = 0.160f;
+ float diffuserLen = 0.250f, diffuserEndDia = 0.070f; // belly
+ float bellyLen = 0.240f;
+ float convergentLen = 0.120f;
+ float stingerDia = 0.026f, stingerLen = 0.070f;
+ // total = 0.16 + 0.25 + 0.24 + 0.12 + 0.07 = 0.84 m
+
+ exhaustHeaderArea = MathF.PI * 0.25f * headerDia * headerDia;
+ float bellyArea = MathF.PI * 0.25f * diffuserEndDia * diffuserEndDia;
+ float stingerArea = MathF.PI * 0.25f * stingerDia * stingerDia;
+
+ float totalExhaustLen = headerLen + diffuserLen + bellyLen + convergentLen + stingerLen; // 840 mm
+ int headerCells = (int)(exhaustCells * (headerLen / totalExhaustLen));
+ int diffuserCells = (int)(exhaustCells * (diffuserLen / totalExhaustLen));
+ int bellyCells = (int)(exhaustCells * (bellyLen / totalExhaustLen));
+ int convergentCells = (int)(exhaustCells * (convergentLen / totalExhaustLen));
+ int stingerCells = exhaustCells - headerCells - diffuserCells - bellyCells - convergentCells;
+
+ // Fill cells
+ for (int i = 0; i < intakeCells; i++)
+ { area[i] = intakePipeArea; dx[i] = intakeLenBefore / intakeCells; }
+ for (int i = intakeCells; i < intakeCells + runnerCells; i++)
+ { area[i] = intakePipeArea; dx[i] = intakeLenRunner / runnerCells; }
+
+ int exhStart = intakeCells + runnerCells;
+ int idx = 0;
+ for (int i = exhStart; i < totalCells; i++)
+ {
+ if (idx < headerCells)
+ { area[i] = exhaustHeaderArea; dx[i] = headerLen / headerCells; }
+ else if (idx < headerCells + diffuserCells)
+ {
+ float t = (idx - headerCells) / (float)(diffuserCells - 1);
+ float dia = headerDia + (diffuserEndDia - headerDia) * t;
+ 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);
+ float dia = diffuserEndDia + (stingerDia - diffuserEndDia) * t;
+ area[i] = MathF.PI * 0.25f * dia * dia;
+ dx[i] = convergentLen / convergentCells;
+ }
+ else
+ { area[i] = stingerArea; dx[i] = stingerLen / stingerCells; }
+ idx++;
+ }
+
+ pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
+ 1.225f, 0f, 101325f);
+ pipeSystem.DampingMultiplier = 1.0f;
+ pipeSystem.EnergyRelaxationRate = 0.5f;
+ pipeSystem.AmbientPressure = 101325f;
+
+ // ---- Volumes ----
+ intakePlenum = new Volume0D(0.5e-3f, 101325f, 300f);
+ plenumInlet = intakePlenum.CreatePort();
+ plenumOutlet = intakePlenum.CreatePort();
+
+ exhaustMuffler = new Volume0D(5e-4f, 101325f, 600f);
+ mufflerIn = exhaustMuffler.CreatePort();
+ mufflerOut = exhaustMuffler.CreatePort();
+
+ // ---- Boundary system ----
+ boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
+ throttleAreaIdx = 0; plenumRunnerIdx = 1; intakeValveIdx = 2; exhaustValveIdx = 3;
+
+ boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, intakePipeArea);
+ intakeOpenIdx = 0;
+ boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, stingerArea);
+ exhaustOpenIdx = 1;
+
+ boundaries.AddOrifice(plenumInlet, 0, false, throttleAreaIdx, 0.7f);
+ boundaries.AddOrifice(plenumOutlet, 1, true, plenumRunnerIdx, 1.0f);
+ boundaries.AddOrifice(cylinder.IntakePort, 1, false, intakeValveIdx, 0.65f);
+ boundaries.AddOrifice(cylinder.ExhaustPort,2, true, exhaustValveIdx, 0.68f);
+
+ orificeAreas = new float[4];
+ orificeAreas[plenumRunnerIdx] = intakePipeArea;
+
+ // ---- Solver ----
+ solver = new Solver { SubStepCount = 4, EnableProfiling = false }; // 4 sub‑steps for 60 cells
+ solver.SetTimeStep(dt);
+ solver.SetPipeSystem(pipeSystem);
+ solver.SetBoundarySystem(boundaries);
+ solver.AddComponent(cylinder);
+ solver.AddComponent(intakePlenum);
+ solver.AddComponent(exhaustMuffler);
+
+ // ---- Sound ----
+ exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
+ intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
+ reverb = new OutdoorExhaustReverb(sampleRate);
+
+ stepCount = 0;
+ Console.WriteLine("125cc Two‑Stroke with vehicle coupling ready.");
+ }
+
+ public override float Process()
+ {
+ float engineRpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+
+ vehicle.ClutchInput = Clutch;
+
+ var (clutchTorque, effectiveInertia) = vehicle.Update(engineRpm, crankshaft.Inertia, (float)dt);
+ crankshaft.SetEffectiveInertia(effectiveInertia);
+ crankshaft.SetLoadTorque(clutchTorque); // clutch torque now includes drag when locked
+
+ 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);
+
+ 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}, Gear={vehicle.CurrentGear}, Speed={vehicle.SpeedKmh:F0} km/h");
+ }
+
+ 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 intakeY = winH / 2f - 40f;
+ float exhaustY = winH / 2f + 80f;
+ float openEndX = 40f;
+
+ // Intake pipe
+ float x = openEndX;
+ float w = 120f;
+ DrawPipe(target, pipeSystem, 0, intakeY, x, x + w);
+
+ // Throttle
+ float throttleX = x + w + 5f;
+ var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
+ {
+ FillColor = Color.Yellow,
+ Position = new Vector2f(throttleX, intakeY - 15f)
+ };
+ target.Draw(throttleRect);
+
+ // Plenum
+ float plenW = 40f, plenH = 60f;
+ float plenX = throttleX + 10f;
+ DrawVolume(target, intakePlenum, plenX + plenW / 2f, intakeY - plenH / 2f, plenW, plenH);
+
+ // Runner
+ float runnerStartX = plenX + plenW + 5f;
+ DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerStartX + 100f);
+
+ // Cylinder
+ float cylCX = runnerStartX + 150f;
+ float cylTopY = intakeY - 120f;
+ DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f);
+
+ // Exhaust pipe
+ float exhStartX = cylCX + 40f + 20f;
+ DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, winW - 60f, areaScale: 1000f);
+
+ // Labels
+ 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);
+ DrawLabel(target, $"Gear: {vehicle.CurrentGear}", new Vector2f(20, 140), Color.Cyan, 20);
+ DrawLabel(target, $"Speed: {vehicle.SpeedKmh:F0} km/h", new Vector2f(20, 160), Color.Cyan, 20);
+
+ // Dyno curve
+ float torqueNm = crankshaft.AverageTorque;
+ UpdateDynoCurve(rpm, powerKw, torqueNm);
+ DrawDynoCurve(target, winW - 410f, winH - 260f, 400f, 250f, rpm, powerKw);
+
+ string gearText = vehicle.CurrentGear == 0 ? "N" : vehicle.CurrentGear.ToString();
+ DrawLabel(target, $"Gear: {gearText}", new Vector2f(20, 140), Color.Cyan, 20);
+ DrawLabel(target, $"Speed: {vehicle.SpeedKmh:F0} km/h", new Vector2f(20, 160), Color.Cyan, 20);
+ DrawLabel(target, vehicle.Engagement > 0.99f ? "Clutch Locked" : "Clutch Slipping", new Vector2f(20, 180), Color.Cyan, 14);
+ }
+ }
+}
\ No newline at end of file