grillkol/FluidSim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: grillkol:a9e1c966b598eb072b24ed52485b05ee2f2269fe
Branches Tags
grillkol:master grillkol:Testing grillkol:testing
..
compare: grillkol:Testing
Branches Tags
grillkol:Testing grillkol:master grillkol:testing

6 Commits
a9e1c966b5 ... Testing

Author SHA1 Message Date
max
56e9c2867a "better" two stroke engine 2026-06-09 22:22:19 +02:00
max
1240ebc33d added two stroke scenario with vehicle 2026-06-09 21:35:48 +02:00
max
ac2eab6f83 250cc mx engine, and dyno 2026-06-09 20:20:56 +02:00
max
aba9b76530 config tuning 2026-06-09 18:05:39 +02:00
max
5c2a7048c8 Merge branch 'Testing' of https://gitea.grillkol.net/grillkol/FluidSim into Testing 2026-06-09 17:50:16 +02:00
max
21a62fb46e stable 2026-06-09 17:49:11 +02:00
12 changed files with 1780 additions and 567 deletions
Expand all files Collapse all files

71
Components/Crankshaft.cs
View File

@@ -4,25 +4,51 @@ 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;
public float FrictionConstant; // N·m
public float FrictionViscous; // N·m per rad/s
public float FrictionConstant;
public float FrictionViscous;
public float LastNetTorque { get; private set; }
public float AveragePower { get; private set; }
public float AverageTorque { get; private set; }
private float externalTorque;
private float _loadTorque;
private readonly float[] _powerBuffer;
private int _powerBufIdx, _powerBufCount;
private float _powerBufSum;
private readonly float[] _torqueBuffer;
private int _torqueBufIdx, _torqueBufCount;
private float _torqueBufSum;
/// <summary>Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke.</summary>
public float CycleLength { get; set; } = 4f * MathF.PI;
public Crankshaft(float initialRPM = 400f)
{
AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
CrankAngle = 0f;
PreviousAngle = 0f;
_powerBuffer = new float[16384];
_torqueBuffer = new float[16384];
}
public void AddTorque(float torque) => externalTorque += torque;
public void SetLoadTorque(float torque) => _loadTorque = Math.Max(torque, 0f);
private float _effectiveInertia; // if >0, overrides Inertia
public void SetEffectiveInertia(float inertia)
{
_effectiveInertia = inertia;
}
public void Step(float dt)
{
if (float.IsNaN(AngularVelocity) || float.IsInfinity(AngularVelocity))
@@ -34,17 +60,42 @@ namespace FluidSim.Components
float friction = FrictionConstant * MathF.Sign(AngularVelocity)
+ FrictionViscous * AngularVelocity;
float netTorque = externalTorque - friction;
float alpha = netTorque / Inertia;
AngularVelocity += alpha * dt;
float netTorque = externalTorque - friction;
LastNetTorque = netTorque;
float totalNetTorque = netTorque - _loadTorque;
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
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
if (_torqueBufCount == _torqueBuffer.Length)
_torqueBufSum -= _torqueBuffer[_torqueBufIdx];
else
_torqueBufCount++;
_torqueBuffer[_torqueBufIdx] = netTorque;
_torqueBufSum += netTorque;
_torqueBufIdx = (_torqueBufIdx + 1) % _torqueBuffer.Length;
AverageTorque = _torqueBufSum / _torqueBufCount;
externalTorque = 0f;
}

210
Components/Cylinder.cs
View File

@@ -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<Port> 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.01f;
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,19 +27,12 @@ namespace FluidSim.Components
float deg = thetaDeg % 720f;
if (deg < 0f) deg += 720f;
float duration;
float effectiveOpen = opens;
float effectiveClose = closes;
if (closes < opens)
{
// Wrap‑around case (e.g., exhaust: opens near 480°, closes near 30°)
effectiveClose += 720f;
}
duration = effectiveClose - effectiveOpen;
if (closes < opens) effectiveClose += 720f;
float duration = effectiveClose - effectiveOpen;
if (duration <= 0f) return 0f;
// Map the angle into the [opens, opens+duration] window
float mapped = deg;
if (mapped < opens) mapped += 720f;
if (mapped < opens || mapped > effectiveClose) return 0f;
@@ -138,39 +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)
{
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
// Intake closing → fuel injection
if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
{
trappedAirMass = _airMass;
@@ -178,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)
@@ -197,7 +89,7 @@ namespace FluidSim.Components
}
}
// Combustion
// Combustion progression
if (combustionActive)
{
float angleSinceSpark = currDeg - sparkAngle;
@@ -220,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;
}
}
}

203
Components/EngineCylinder.cs Normal file
View File

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

183
Components/TwoStrokeCylinder.cs Normal file
View File

@@ -0,0 +1,183 @@
using System;
namespace FluidSim.Components
{
/// <summary>
/// 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.
/// </summary>
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
// ── Configurable durations ──────────────────────────────────────────────
public float TransferDuration { get; } // default: 155°
public float ExhaustDuration { get; } // default: 195°
// 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;
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);
// ── 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;
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.");
}
// ── Valve lift profile ──────────────────────────────────────────────────
/// <summary>
/// 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).
/// </summary>
private static 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;
if (duration <= 0f) return 0f;
// Map deg into the same number-line as opens/effectiveClose
float mapped = deg < opens ? deg + 360f : deg;
if (mapped < opens || mapped > effectiveClose) return 0f;
float rampDur = duration * RampFraction;
float holdEnd = effectiveClose - rampDur;
if (mapped < opens + rampDur)
{
// Opening ramp: smoothstep
float t = (mapped - opens) / rampDur;
return peakLift * t * t * (3f - 2f * t);
}
else if (mapped <= holdEnd)
{
// Flat top – full area
return peakLift;
}
else
{
// Closing ramp: smoothstep reversed
float t = (mapped - holdEnd) / rampDur;
return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
}
}
// ── 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))
{
trappedAirMass = _airMass;
fuelMass = trappedAirMass / StoichiometricAFR;
fuelInjected = true;
}
// ── Ignition ───────────────────────────────────────────────────────
// SparkAdvance default is ~22° BTDC on the base class; scenario can override.
float sparkAngle = (360f - SparkAdvance) % 360f;
if (CrossedAngle(prevDeg, currDeg, sparkAngle) && !combustionActive && fuelInjected)
{
if (_random.NextDouble() < MisfireProbability)
{
combustionActive = false;
}
else
{
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)
{
newFraction = 1f;
combustionActive = false;
fuelInjected = false;
float totalMass = _airMass + _exhaustMass;
_airMass = 0f;
_exhaustMass = totalMass;
}
float dFraction = newFraction - burnFraction;
if (dFraction > 0f)
{
float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction;
cylinderEnergy += dQ;
_exhaustMass += fuelMass * dFraction;
burnFraction = newFraction;
}
}
}
// ── Helper: did the crank cross a target angle this step? ───────────────
/// <summary>
/// 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)
if (curr >= prev)
return prev < target && target <= curr;
// Wrapped past 360° → two intervals to check
return prev < target || target <= curr;
}
}
}

166
Components/Vehicle.cs Normal file
View File

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

144
Core/BoundarySystem.cs
View File

@@ -25,8 +25,7 @@ namespace FluidSim.Core
public float EffectiveLength;
public float CurrentMdot; // kg/s, positive = volume → pipe
// --- Loss coefficient (linear resistance) – inertance only ---
// If 0 when UseInertance is true, a stable default is auto‑computed at runtime.
// --- Loss coefficient (linear resistance) ---
public float LossCoefficient; // N·s/m⁵ or kg/(m⁴·s)
}
@@ -58,10 +57,9 @@ namespace FluidSim.Core
public int OpenEndCount { get; private set; }
// ---------- Add orifice (no inertance) ----------
// Simple isentropic nozzle – no built‑in loss. For dissipation use pipe damping
// or the inertance model if you need a damped resonator.
public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
int areaIndex, float dischargeCoeff = 1f)
int areaIndex, float dischargeCoeff = 1f,
float lossCoefficient = 0f)
{
_orifices[OrificeCount] = new OrificeDesc
{
@@ -73,24 +71,22 @@ namespace FluidSim.Core
UseInertance = false,
EffectiveLength = 0f,
CurrentMdot = 0f,
LossCoefficient = 0f
LossCoefficient = lossCoefficient
};
OrificeCount++;
}
// ---------- Add orifice with inertance ----------
// effectiveLength – length of the inertial slug (m), typically the physical neck length.
// lossCoefficient – linear resistance (N·s/m⁵). If 0 (or omitted) an automatic stable
// value will be computed from the pipe's characteristic impedance.
public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
int areaIndex, float dischargeCoeff,
float effectiveLength, float lossCoefficient = 0f)
{
AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff);
// Reuse the base AddOrifice and then override fields
AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient);
ref var d = ref _orifices[OrificeCount - 1];
d.UseInertance = true;
d.EffectiveLength = effectiveLength;
d.LossCoefficient = lossCoefficient;
d.LossCoefficient = lossCoefficient; // store the linear resistance
}
public void AddOpenEnd(int pipeIndex, bool isLeftEnd,
@@ -150,7 +146,7 @@ namespace FluidSim.Core
? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
: _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
// ---- Handle closed orifice as a wall ----
// ---- Handle closed orifice (area ≈ 0) as a wall ----
if (area < 1e-12f || d.VolumePort == null)
{
var (rInt, uInt, pInt) = d.IsLeftEnd
@@ -169,7 +165,7 @@ namespace FluidSim.Core
continue;
}
// ---- Preliminary isentropic solution (for reference) ----
// ---- Preliminary isentropic solution ----
float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
if (volP >= pipeP)
{
@@ -183,31 +179,20 @@ namespace FluidSim.Core
mdotEst = -mdotEst;
}
// ---- Compute ghost state ----
// ---- Compute final mass flow with limiters ----
float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
if (d.UseInertance)
{
// ---- Inertance ODE with (possibly automatic) linear loss ----
float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
float dp = volP - pipeP;
// If loss coefficient not provided, use a tiny fraction of the pipe's characteristic impedance
float Rlin = d.LossCoefficient;
if (Rlin <= 0f)
{
// Auto‑sized linear drag: 0.5 % of Z_char
float rhoRef = d.CurrentMdot >= 0 ? volRho : pipeRho;
float cRef = d.CurrentMdot >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
float Z_char = rhoRef * cRef / MathF.Max(area, 1e-12f);
Rlin = 0.005f * Z_char;
}
float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
float mdotNew = d.CurrentMdot + dmdot_dt * dt;
// Symmetric flow limiters
// Limit outflow from volume (if volume owner is Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
@@ -215,15 +200,19 @@ namespace FluidSim.Core
if (mdotNew < -maxOut) mdotNew = -maxOut;
}
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
// Limit inflow from pipe – pipe cell cannot be emptied in one step
{
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); // true cell area, not orifice
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
}
// Velocity clamp Mach 0.9
// Velocity clamp to Mach 0.9
float rhoFacePrelim = mdotNew >= 0 ? volRho : pipeRho;
float uFacePrelim = MathF.Abs(mdotNew) / MathF.Max(rhoFacePrelim * area, 1e-12f);
float cUp = mdotNew >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
@@ -238,51 +227,60 @@ namespace FluidSim.Core
d.CurrentMdot = mdotNew;
mdotFinal = mdotNew;
rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
pFace = pFaceEst;
uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
}
else
{
// ---- Standard quasi‑steady orifice (purely isentropic) ----
// Standard quasi‑steady orifice
mdotFinal = mdotEst;
rhoFace = rhoFaceEst;
uFace = uFaceEst;
pFace = pFaceEst;
// Limit outflow from cavity
// Limit outflow from volume (if Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
if (mdotFinal > maxOut) mdotFinal = maxOut;
}
// Safety velocity clamp (Mach 0.9)
float cLocal = mdotFinal >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
float maxULocal = 0.9f * cLocal;
float uCheck = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
if (uCheck > maxULocal)
// ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
if (mdotFinal < 0)
{
uFace = maxULocal;
mdotFinal = rhoFace * uFace * area * (mdotFinal >= 0 ? 1f : -1f);
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotFinal < -maxFromPipe)
mdotFinal = -maxFromPipe;
}
d.CurrentMdot = mdotFinal;
// Limit outflow from cylinder into pipe (positive mdot = volume → pipe)
if (mdotFinal > 0f && d.VolumePort?.Owner is Cylinder cyl)
{
float maxOut = cyl.Mass / dt;
if (mdotFinal > maxOut)
mdotFinal = maxOut;
}
}
// ---- Determine air fraction for ghost ----
// ---- Air fraction for ghost ----
if (mdotFinal >= 0)
{
airFracGhost = volAF;
}
else
{
airFracGhost = pipeAF;
if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
}
// ---- Apply sign convention for velocity ----
// ---- Sign convention for velocity ----
if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace;
else if (mdotFinal >= 0 && !d.IsLeftEnd) uFace = -uFace;
else if (mdotFinal < 0 && d.IsLeftEnd) uFace = -uFace;
@@ -299,12 +297,12 @@ namespace FluidSim.Core
{
d.VolumePort.MassFlowRate = -mdotFinal;
if (-mdotFinal >= 0) // mass flowing into the volume
if (-mdotFinal >= 0) // mass entering volume (out of pipe)
{
float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
d.VolumePort.SpecificEnthalpy = pipeH;
}
else // mass flowing out of the volume
else // mass leaving volume (into pipe)
{
d.VolumePort.SpecificEnthalpy = volH;
}
@@ -331,6 +329,7 @@ namespace FluidSim.Core
float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f));
float pAmb = d.AmbientPressure;
// Characteristic solution (isentropic expansion to ambient)
float Jplus = uInt + 2f * cInt / gm1;
float Jminus = uInt - 2f * cInt / gm1;
float s = pInt / MathF.Pow(rhoInt, gamma);
@@ -340,9 +339,14 @@ namespace FluidSim.Core
? (Jminus + 2f * cIso / gm1)
: (Jplus - 2f * cIso / gm1);
// Supersonic check
bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
float rhoGhost, uGhost, pGhost, afGhost;
if (!supersonic)
{
supersonic = d.IsLeftEnd ? (uIso <= -cIso) : (uIso >= cIso);
}
float rhoGhost, uGhost, pGhost, afGhost;
if (supersonic)
{
rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt;
@@ -354,15 +358,45 @@ namespace FluidSim.Core
afGhost = inflow ? 1f : afInt;
}
// ------- Mass flow limiter -------
int adjCell = d.IsLeftEnd
? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float cellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float area = d.PipeArea;
float mdotRaw = rhoGhost * uGhost * area; // positive out of pipe
if (d.IsLeftEnd) mdotRaw = -mdotRaw; // now positive = out of pipe
// Outflow limit
if (mdotRaw > 0 && mdotRaw * dt > cellMass)
{
mdotRaw = cellMass / dt;
}
// Inflow limit (allow up to 10× cell mass to avoid starving the pipe)
else if (mdotRaw < 0 && -mdotRaw * dt > 10f * cellMass)
{
mdotRaw = -10f * cellMass / dt;
}
// Recompute uGhost from the limited mdot, keeping rhoGhost, pGhost
float mdotMag = MathF.Abs(mdotRaw);
uGhost = mdotMag / MathF.Max(rhoGhost * area, 1e-12f);
if (d.IsLeftEnd)
uGhost = (mdotRaw >= 0f) ? -uGhost : uGhost;
else
uGhost = (mdotRaw >= 0f) ? uGhost : -uGhost;
// Apply ghost
if (d.IsLeftEnd)
_pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
else
_pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
float area = d.PipeArea;
float mdot = rhoGhost * uGhost * area;
if (d.IsLeftEnd) mdot = -mdot;
d.LastMassFlowRate = mdot;
d.LastMassFlowRate = mdotRaw;
d.LastFacePressure = pGhost;
}
}

428
Core/Pipesystem.cs
View File

@@ -16,23 +16,28 @@ namespace FluidSim.Core
private readonly int _allCells; // total allocated (padded to Vector<float>.Count)
private readonly int _pipeCount;
// Derived state – _p is kept for visualization, _c is gone
// Derived state – _p is kept for visualization
private float[] _p;
// Flux arrays (size = _allCells + 1)
// Flux arrays for faces INTERNAL to a single pipe (size = _allCells + 1)
// Only valid for faces that are NOT pipe boundaries.
private float[] _fluxM, _fluxP, _fluxE, _fluxY;
// Damping and relaxation (computed on‑the‑fly only if used)
// Per‑pipe boundary flux buffers (size = _pipeCount)
private float[] _leftFluxM, _leftFluxP, _leftFluxE, _leftFluxY;
private float[] _rightFluxM, _rightFluxP, _rightFluxE, _rightFluxY;
// Damping and relaxation
private float[] _dampingFactors;
private float[] _relaxFactors;
private bool _applyDamping;
private bool _applyRelax;
// Ghost buffer
// Ghost buffer (per‑pipe ghost states)
private readonly GhostBuffer _ghost;
// Wall mask – precomputed once
private readonly bool[] _isWallFace;
// Precomputed flag: true if a face is a pipe boundary (start or end)
private readonly bool[] _isPipeBoundaryFace;
// ---------- Physical constants ----------
private const float Gamma = 1.4f;
@@ -102,6 +107,16 @@ namespace FluidSim.Core
_fluxE = new float[faceCount];
_fluxY = new float[faceCount];
// Per‑pipe boundary flux buffers
_leftFluxM = new float[_pipeCount];
_leftFluxP = new float[_pipeCount];
_leftFluxE = new float[_pipeCount];
_leftFluxY = new float[_pipeCount];
_rightFluxM = new float[_pipeCount];
_rightFluxP = new float[_pipeCount];
_rightFluxE = new float[_pipeCount];
_rightFluxY = new float[_pipeCount];
_dampingFactors = new float[_allCells];
_relaxFactors = new float[_allCells];
_applyDamping = _coeffBase != 0f;
@@ -110,18 +125,12 @@ namespace FluidSim.Core
_ghost = new GhostBuffer(_pipeCount);
_ambientEnergyRef = initialP * Gm1Inv;
// Pre‑compute wall face flags: each face that sits between two different pipes is a wall
_isWallFace = new bool[faceCount];
for (int f = 1; f < _totalCells; f++)
// Mark faces that coincide with a pipe boundary (start or end)
_isPipeBoundaryFace = new bool[faceCount];
for (int p = 0; p < _pipeCount; p++)
{
for (int p = 0; p < _pipeCount; p++)
{
if (f == _pipeEnd[p] && f < _totalCells)
{
_isWallFace[f] = true;
break;
}
}
_isPipeBoundaryFace[_pipeStart[p]] = true;
_isPipeBoundaryFace[_pipeEnd[p]] = true;
}
// Initialize uniform state
@@ -150,6 +159,7 @@ namespace FluidSim.Core
public float GetCellPressure(int i) => _p[i];
public float GetCellDensity(int i) => _rho[i];
public float GetCellDx(int i) => _dx[i];
public float GetCellArea(int i) => _area[i];
public float GetCellVelocity(int i)
{
float rho = _rho[i];
@@ -215,13 +225,13 @@ namespace FluidSim.Core
}
}
// ---------- Flux computation: fuses primitive calculation and flux evaluation ----------
// ---------- Flux computation ----------
private void ComputeFluxes(float dt)
{
float fm, fp, fe;
int vecSize = Vector<float>.Count;
// ---- 1. Left ghost boundaries ----
// ---- 1. Left ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeStart[p];
@@ -239,22 +249,18 @@ namespace FluidSim.Core
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
float YR = _Y[idx];
// store pressure for cell idx
_p[idx] = pR;
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[idx] = fm; _fluxP[idx] = fp; _fluxE[idx] = fe;
_leftFluxM[p] = fm; _leftFluxP[p] = fp; _leftFluxE[p] = fe;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_fluxY[idx] = fy;
_leftFluxY[p] = fy;
}
// ---- 2. Right ghost boundaries ----
// ---- 2. Right ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeEnd[p] - 1;
int face = idx + 1;
int ghostIdx = p * 2 + 1;
float rR = _ghost.Rho[ghostIdx];
float uR = _ghost.U[ghostIdx];
@@ -269,45 +275,35 @@ namespace FluidSim.Core
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
float YL = _Y[idx];
// store pressure for cell idx
_p[idx] = pL;
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
_rightFluxM[p] = fm; _rightFluxP[p] = fp; _rightFluxE[p] = fe;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_fluxY[face] = fy;
_rightFluxY[p] = fy;
}
// ---- 3. Interior faces – vectorised SIMD ----
// ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
for (int face = 1; face < _totalCells; face++)
{
// Handle walls (rare) with scalar code
if (_isWallFace[face])
{
int iL = face - 1;
float rL = _rho[iL], rhouL = _rhou[iL];
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
float uL = rhouL * invRhoL;
float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
_p[iL] = pL;
LaxFlux(rL, uL, pL, cL, rL, -uL, pL, cL, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
_fluxY[face] = 0f;
// Skip faces that belong to a pipe boundary (they are already handled)
if (_isPipeBoundaryFace[face])
continue;
}
// If the next vecSize faces contain a wall, fall back to scalar for this block
// Try to vectorize a block of contiguous non‑boundary faces
if (face + vecSize - 1 < _totalCells)
{
bool hasWall = false;
bool canVectorize = true;
for (int f = face; f < face + vecSize; f++)
if (_isWallFace[f]) { hasWall = true; break; }
{
if (_isPipeBoundaryFace[f])
{
canVectorize = false;
break;
}
}
if (!hasWall)
if (canVectorize)
{
// --- Vectorised block ---
var rhoL = new Vector<float>(_rho, face - 1);
@@ -330,11 +326,7 @@ namespace FluidSim.Core
var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
// Store pressures for visualisation (left cell of each face)
pL.CopyTo(_p, face - 1);
// Lax‑Friedrichs fluxes
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
var FmL = rhoL * uL;
@@ -362,50 +354,45 @@ namespace FluidSim.Core
feVec.CopyTo(_fluxE, face);
fyVec.CopyTo(_fluxY, face);
face += vecSize - 1; // loop increment will add 1, so we advance vecSize faces
face += vecSize - 1; // loop increment will add 1
continue;
}
}
// --- Scalar interior face (fallback) ---
// --- Scalar fallback for a single interior face ---
{
int iLf = face - 1, iRf = face;
float rLf = _rho[iLf], rhouLf = _rhou[iLf];
float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f));
float uLf = rhouLf * invRhoLf;
float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf);
float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf);
float YLf = _Y[iLf];
_p[iLf] = pLf;
int iL = face - 1, iR = face;
float rL = _rho[iL], rhouL = _rhou[iL];
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
float uL = rhouL * invRhoL;
float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
float YL = _Y[iL];
float rRf = _rho[iRf], rhouRf = _rhou[iRf];
float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f));
float uRf = rhouRf * invRhoRf;
float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf);
float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf);
float YRf = _Y[iRf];
float rR = _rho[iR], rhouR = _rhou[iR];
float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
float uR = rhouR * invRhoR;
float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
float YR = _Y[iR];
LaxFlux(rLf, uLf, pLf, cLf, rRf, uRf, pRf, cRf, out fm, out fp, out fe);
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
float alpha = MathF.Max(MathF.Abs(uLf) + cLf, MathF.Abs(uRf) + cRf);
ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy);
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_fluxY[face] = fy;
}
}
// If damping/relaxation are active, compute the factors here (re-uses _dampingFactors/_relaxFactors arrays,
// but we no longer have a separate precompute pass). We compute them on demand in UpdateCells anyway?
// Actually UpdateCells multiplies by these factors; we can compute them there if needed.
}
// ---------- Cell update (unchanged core, but skips relaxation/damping when not needed) ----------
// ---------- Cell update (per pipe, using correct boundary fluxes) ----------
private void UpdateCells(float dt)
{
int vecSize = Vector<float>.Count;
float dtRelax = -_relaxRate * dt;
// Compute damping and relaxation factors if needed
// Precompute damping and relaxation factors globally
if (_applyDamping)
{
for (int i = 0; i < _totalCells; i++)
@@ -418,89 +405,217 @@ namespace FluidSim.Core
}
if (_applyRelax)
{
var relaxVal = MathF.Exp(dtRelax);
float relaxVal = MathF.Exp(dtRelax);
for (int i = 0; i < _totalCells; i++)
_relaxFactors[i] = relaxVal;
}
int iCell = 0;
for (; iCell <= _totalCells - vecSize; iCell += vecSize)
// Update each pipe separately
for (int p = 0; p < _pipeCount; p++)
{
var rhoOld = new Vector<float>(_rho, iCell);
var rhouOld = new Vector<float>(_rhou, iCell);
var EOld = new Vector<float>(_E, iCell);
var YOld = new Vector<float>(_Y, iCell);
int start = _pipeStart[p];
int end = _pipeEnd[p]; // exclusive
int len = end - start;
if (len == 0) continue;
var fluxM_L = new Vector<float>(_fluxM, iCell);
var fluxP_L = new Vector<float>(_fluxP, iCell);
var fluxE_L = new Vector<float>(_fluxE, iCell);
var fluxY_L = new Vector<float>(_fluxY, iCell);
var fluxM_R = new Vector<float>(_fluxM, iCell + 1);
var fluxP_R = new Vector<float>(_fluxP, iCell + 1);
var fluxE_R = new Vector<float>(_fluxE, iCell + 1);
var fluxY_R = new Vector<float>(_fluxY, iCell + 1);
var dtdx = new Vector<float>(dt) / new Vector<float>(_dx, iCell);
var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
var ENew = EOld - dtdx * (fluxE_R - fluxE_L);
var rhoYOld = rhoOld * YOld;
var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping)
rhouNew *= new Vector<float>(_dampingFactors, iCell);
if (_applyRelax)
// ------- Left boundary cell (i = start) ------
{
var ambRef = new Vector<float>(_ambientEnergyRef);
var relax = new Vector<float>(_relaxFactors, iCell);
ENew = ambRef + (ENew - ambRef) * relax;
int i = start;
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
// left face: always the pipe's left boundary flux
float fluxM_L = _leftFluxM[p];
float fluxP_L = _leftFluxP[p];
float fluxE_L = _leftFluxE[p];
float fluxY_L = _leftFluxY[p];
// right face: depends on pipe length
float fluxM_R, fluxP_R, fluxE_R, fluxY_R;
if (len == 1)
{
// Only one cell: right face is the pipe's right boundary flux
fluxM_R = _rightFluxM[p];
fluxP_R = _rightFluxP[p];
fluxE_R = _rightFluxE[p];
fluxY_R = _rightFluxY[p];
}
else
{
// interior face (global flux at index i+1)
fluxM_R = _fluxM[i + 1];
fluxP_R = _fluxP[i + 1];
fluxE_R = _fluxE[i + 1];
fluxY_R = _fluxY[i + 1];
}
float dtdx = dt / _dx[i];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[i];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[i] = rhoNew;
_rhou[i] = rhouNew;
_E[i] = ENew;
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
rhoNew = Vector.Max(rhoNew, new Vector<float>(1e-12f));
var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
ENew = Vector.Max(ENew, eMin);
// ------- Interior cells (i = start+1 to end-2) ------
if (len > 2)
{
int iCell = start + 1;
int iEnd = end - 1; // exclusive upper bound
rhoNew.CopyTo(_rho, iCell);
rhouNew.CopyTo(_rhou, iCell);
ENew.CopyTo(_E, iCell);
var yNew = rhoYNew / rhoNew;
yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
yNew.CopyTo(_Y, iCell);
// Vectorised path for interior cells (if available)
for (; iCell <= iEnd - vecSize; iCell += vecSize)
{
var rhoOld = new Vector<float>(_rho, iCell);
var rhouOld = new Vector<float>(_rhou, iCell);
var EOld = new Vector<float>(_E, iCell);
var YOld = new Vector<float>(_Y, iCell);
var fluxM_L = new Vector<float>(_fluxM, iCell);
var fluxP_L = new Vector<float>(_fluxP, iCell);
var fluxE_L = new Vector<float>(_fluxE, iCell);
var fluxY_L = new Vector<float>(_fluxY, iCell);
var fluxM_R = new Vector<float>(_fluxM, iCell + 1);
var fluxP_R = new Vector<float>(_fluxP, iCell + 1);
var fluxE_R = new Vector<float>(_fluxE, iCell + 1);
var fluxY_R = new Vector<float>(_fluxY, iCell + 1);
var dtdx = new Vector<float>(dt) / new Vector<float>(_dx, iCell);
var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
var ENew = EOld - dtdx * (fluxE_R - fluxE_L);
var rhoYOld = rhoOld * YOld;
var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping)
rhouNew *= new Vector<float>(_dampingFactors, iCell);
if (_applyRelax)
{
var ambRef = new Vector<float>(_ambientEnergyRef);
var relax = new Vector<float>(_relaxFactors, iCell);
ENew = ambRef + (ENew - ambRef) * relax;
}
rhoNew = Vector.Max(rhoNew, new Vector<float>(1e-12f));
var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
ENew = Vector.Max(ENew, eMin);
rhoNew.CopyTo(_rho, iCell);
rhouNew.CopyTo(_rhou, iCell);
ENew.CopyTo(_E, iCell);
var yNew = rhoYNew / rhoNew;
yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
yNew.CopyTo(_Y, iCell);
}
// Scalar remainder for interior cells
for (; iCell < iEnd; iCell++)
{
float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
float dtdx = dt / _dx[iCell];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[iCell];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[iCell] = rhoNew;
_rhou[iCell] = rhouNew;
_E[iCell] = ENew;
_Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
}
// ------- Right boundary cell (i = end-1, if len > 1) ------
if (len > 1)
{
int i = end - 1;
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
// left face
float fluxM_L, fluxP_L, fluxE_L, fluxY_L;
if (len == 2)
{
// Only two cells: left face is the pipe's left boundary flux
fluxM_L = _leftFluxM[p];
fluxP_L = _leftFluxP[p];
fluxE_L = _leftFluxE[p];
fluxY_L = _leftFluxY[p];
}
else
{
// interior face (global flux at i)
fluxM_L = _fluxM[i];
fluxP_L = _fluxP[i];
fluxE_L = _fluxE[i];
fluxY_L = _fluxY[i];
}
// right face: always the pipe's right boundary flux
float fluxM_R = _rightFluxM[p];
float fluxP_R = _rightFluxP[p];
float fluxE_R = _rightFluxE[p];
float fluxY_R = _rightFluxY[p];
float dtdx = dt / _dx[i];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[i];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[i] = rhoNew;
_rhou[i] = rhouNew;
_E[i] = ENew;
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
}
// Scalar remainder (only a few cells)
for (; iCell < _totalCells; iCell++)
// Recompute pressure for all cells (for visualization)
for (int i = 0; i < _totalCells; i++)
{
float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
float dtdx = dt / _dx[iCell];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[iCell];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[iCell] = rhoNew;
_rhou[iCell] = rhouNew;
_E[iCell] = ENew;
_Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
float rho = _rho[i];
float rhou = _rhou[i];
float u = rhou / MathF.Max(rho, 1e-12f);
_p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
}
}
// ---------- Scalar flux helpers (used in boundaries and scalar fallback) ----------
// ---------- Scalar flux helpers ----------
private static void LaxFlux(float rL, float uL, float pL, float cL,
float rR, float uR, float pR, float cR,
out float fm, out float fp, out float fe)
@@ -528,6 +643,23 @@ namespace FluidSim.Core
fy = 0.5f * (FyL + FyR) - 0.5f * alpha * (rR * YR - rL * YL);
}
public int GetRequiredSubSteps(float dtGlobal, float cflTarget = 0.8f)
{
float maxW = 0f;
for (int i = 0; i < _totalCells; i++)
{
float rho = MathF.Max(_rho[i], 1e-12f);
float u = MathF.Abs(_rhou[i] / rho);
float p = Gm1 * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
float c = MathF.Sqrt(Gamma * p / rho);
float w = u + c;
if (w > maxW) maxW = w;
}
maxW = MathF.Max(maxW, 1e-8f);
float minDx = _dx.Min(); // need using System.Linq;
return Math.Max(1, (int)MathF.Ceiling(dtGlobal * maxW / (cflTarget * minDx)));
}
// ---------- Profiling report ----------
public string GetProfileReport()
{

3
Core/Solver.cs
View File

@@ -36,7 +36,8 @@ namespace FluidSim.Core
{
if (_pipeSystem == null || _boundarySystem == null) return;
int nSub = SubStepCount;
int nSub = _pipeSystem.GetRequiredSubSteps((float)_dt, 0.8f);
nSub = Math.Max(nSub, SubStepCount); // never go below fixed minimum
float dtSub = (float)(_dt / nSub);
for (int sub = 0; sub < nSub; sub++)

60
Program.cs
View File

@@ -33,24 +33,33 @@ public class Program
// Audio & simulation
private static SimulationRingBuffer _simRingBuffer = null!;
private static SoundEngine _soundEngine = null!;
private static Scenario _scenario = null!; // cast to access ThrottleArea
private static Scenario _scenario = null!;
private static Font? _overlayFont;
private static Text? _overlayText;
// Throttle control
private static float _throttleTarget = 1.0f; // 0‑1, set by arrow keys
private static float _throttleCurrent = 0.0f; // actual current fraction (lerped)
private const float ThrottleLerpRate = 10.0f; // times per second (speed of movement)
private static float _throttleTarget = 1.0f;
private static float _throttleCurrent = 0.0f;
private const float ThrottleLerpRate = 10.0f;
private static bool _wKeyHeld = false;
private static float _lastThrottleUpdateTime;
// Load
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;
@@ -76,14 +85,12 @@ public class Program
(1.0 - Math.Exp(-8.0 * (now - lastDrawTime)));
_soundEngine.Speed = _currentDisplaySpeed;
// ---- Throttle update ----
// ---- Throttle & Load update (shared dt) ----
float dtThrottle = (float)now - _lastThrottleUpdateTime;
_lastThrottleUpdateTime = (float)now;
float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
// Snap to zero instantly when target is zero (key released)
if (throttleDesiredFraction == 0.0)
if (throttleDesiredFraction == 0.0f)
{
_throttleCurrent = 0.0f;
}
@@ -93,8 +100,18 @@ public class Program
_throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
}
float loadSmoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
_loadCurrent += (_loadTarget - _loadCurrent) * loadSmoothing;
_scenario.Load = _loadCurrent;
_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)
{
@@ -103,7 +120,8 @@ 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" +
$"Throttle: {_throttleCurrent * 100:F0}% Target: {_throttleTarget * 100:F0}% [W] {(_wKeyHeld ? "BLIP" : "---")}";
$"Clutch: {_clutchCurrent*100:F0}% [C]" +
$"Load: {_loadCurrent*100:F0}% [←][→] Throttle: {_throttleCurrent * 100:F0}% Target: {_throttleTarget * 100:F0}% [W] {(_wKeyHeld ? "BLIP" : "---")}";
}
window.Clear(Color.Black);
@@ -205,6 +223,25 @@ public class Program
case Keyboard.Key.Down:
_throttleTarget = MathF.Max(0.0f, _throttleTarget - 0.05f);
break;
case Keyboard.Key.Left:
_loadTarget = MathF.Max(0.0f, _loadTarget - 0.05f);
break;
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;
}
}
@@ -212,5 +249,8 @@ public class Program
{
if (e.Code == Keyboard.Key.W)
_wKeyHeld = false;
if (e.Code == Keyboard.Key.C)
_cKeyHeld = false;
}
}

235
Scenarios/Scenario.cs
View File

@@ -2,6 +2,8 @@
using SFML.System;
using FluidSim.Core;
using FluidSim.Components;
using System;
using System.Collections.Generic;
namespace FluidSim.Tests
{
@@ -10,11 +12,204 @@ namespace FluidSim.Tests
protected const float AmbientPressure = 101325f;
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();
private int _lastDynoBin = -1;
public void ResetDynoCurve()
{
_dynoBins.Clear();
_lastDynoBin = -1;
}
protected void UpdateDynoCurve(float rpm, float powerKw, float torqueNm)
{
if (rpm <= 0) return;
int bin = (int)(rpm / RpmBinSize);
while (_dynoBins.Count <= bin)
_dynoBins.Add((0f, 0f));
if (_lastDynoBin >= 0 && bin > _lastDynoBin + 1)
{
float lastPower = _dynoBins[_lastDynoBin].powerKw > 0 ? _dynoBins[_lastDynoBin].powerKw : 0f;
float lastTorque = _dynoBins[_lastDynoBin].torqueNm > 0 ? _dynoBins[_lastDynoBin].torqueNm : 0f;
for (int b = _lastDynoBin + 1; b < bin; b++)
{
float t = (b - _lastDynoBin) / (float)(bin - _lastDynoBin);
float interpPower = lastPower + (powerKw - lastPower) * t;
float interpTorque = lastTorque + (torqueNm - lastTorque) * t;
if (interpPower > _dynoBins[b].powerKw || _dynoBins[b].powerKw <= 0)
_dynoBins[b] = (interpPower, _dynoBins[b].torqueNm);
if (interpTorque > _dynoBins[b].torqueNm || _dynoBins[b].torqueNm <= 0)
_dynoBins[b] = (_dynoBins[b].powerKw, interpTorque);
}
}
var current = _dynoBins[bin];
if (powerKw > current.powerKw || current.powerKw <= 0)
current.powerKw = powerKw;
if (torqueNm > current.torqueNm || current.torqueNm <= 0)
current.torqueNm = torqueNm;
_dynoBins[bin] = current;
_lastDynoBin = bin;
}
protected void DrawDynoCurve(RenderWindow target,
float graphX, float graphY, float graphWidth, float graphHeight,
float currentRpm, float currentPowerKw)
{
if (_dynoBins.Count == 0) return;
float maxPowerKw = 0.01f, maxTorqueNm = 0.01f, maxRpm = 1000f;
for (int b = 0; b < _dynoBins.Count; b++)
{
var bin = _dynoBins[b];
if (bin.powerKw > 0 || bin.torqueNm > 0)
{
float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
if (bin.powerKw > maxPowerKw) maxPowerKw = bin.powerKw;
if (bin.torqueNm > maxTorqueNm) maxTorqueNm = bin.torqueNm;
if (rpmBin > maxRpm) maxRpm = rpmBin;
}
}
maxPowerKw *= 1.1f;
maxTorqueNm *= 1.1f;
maxRpm = MathF.Max(maxRpm * 1.05f, 1000f);
var bg = new RectangleShape(new Vector2f(graphWidth, graphHeight))
{
FillColor = new Color(20, 20, 20, 200),
Position = new Vector2f(graphX, graphY)
};
target.Draw(bg);
const float leftMargin = 50f, rightMargin = 50f, topMargin = 20f, bottomMargin = 35f;
float plotX = graphX + leftMargin;
float plotY = graphY + topMargin;
float plotW = graphWidth - leftMargin - rightMargin;
float plotH = graphHeight - topMargin - bottomMargin;
float xMin = 0f, xMax = maxRpm;
float yLeftMin = 0f, yLeftMax = maxPowerKw;
float yRightMin = 0f, yRightMax = maxTorqueNm;
var powerColor = new Color(0xFF, 0x1B, 0x1B);
var torqueColor = new Color(0x09, 0x09, 0xFF);
var gridColor = new Color(50, 50, 50);
for (int i = 0; i <= 9; i++)
{
float t = i / 9f;
float x = plotX + t * plotW;
var vLine = new VertexArray(PrimitiveType.Lines, 2);
vLine[0] = new Vertex(new Vector2f(x, plotY), gridColor);
vLine[1] = new Vertex(new Vector2f(x, plotY + plotH), gridColor);
target.Draw(vLine);
}
for (int i = 0; i <= 5; i++)
{
float t = i / 5f;
float y = plotY + (1 - t) * plotH;
var hLine = new VertexArray(PrimitiveType.Lines, 2);
hLine[0] = new Vertex(new Vector2f(plotX, y), gridColor);
hLine[1] = new Vertex(new Vector2f(plotX + plotW, y), gridColor);
target.Draw(hLine);
}
DrawLabel(target, "RPM", new Vector2f(graphX + graphWidth / 2 - 12, graphY + graphHeight - 15), Color.White, 12);
DrawLabel(target, "kW", new Vector2f(graphX + 5, graphY + 2), Color.White, 11);
DrawLabel(target, "Nm", new Vector2f(graphX + graphWidth - 25, graphY + 2), Color.White, 11);
for (int i = 0; i <= 5; i++)
{
float leftValue = yLeftMin + (yLeftMax - yLeftMin) * i / 5f;
float rightValue = yRightMin + (yRightMax - yRightMin) * i / 5f;
float y = plotY + (1 - i / 5f) * plotH;
DrawLabel(target, $"{leftValue:F1}", new Vector2f(graphX + 2, y - 6), Color.White, 9);
DrawLabel(target, $"{rightValue:F1}", new Vector2f(graphX + graphWidth - 40, y - 6), Color.White, 9);
}
for (int i = 0; i <= 9; i++)
{
float value = xMin + (xMax - xMin) * i / 9f;
float x = plotX + i / 9f * plotW;
DrawLabel(target, $"{value / 1000f:F1}k", new Vector2f(x - 15, graphY + graphHeight - bottomMargin + 5), Color.White, 9);
}
var powerLine = new VertexArray(PrimitiveType.LineStrip);
bool firstPower = true;
for (int b = 0; b < _dynoBins.Count; b++)
{
float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
if (rpmBin > xMax) break;
var bin = _dynoBins[b];
if (bin.powerKw > 0)
{
float sx = plotX + (rpmBin - xMin) / (xMax - xMin) * plotW;
float sy = plotY + (1 - (bin.powerKw - yLeftMin) / (yLeftMax - yLeftMin)) * plotH;
if (firstPower) { powerLine.Clear(); firstPower = false; }
powerLine.Append(new Vertex(new Vector2f(sx, sy), powerColor));
}
else if (!firstPower)
{
target.Draw(powerLine);
powerLine.Clear();
firstPower = true;
}
}
if (!firstPower) target.Draw(powerLine);
var torqueLine = new VertexArray(PrimitiveType.LineStrip);
bool firstTorque = true;
for (int b = 0; b < _dynoBins.Count; b++)
{
float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
if (rpmBin > xMax) break;
var bin = _dynoBins[b];
if (bin.torqueNm > 0)
{
float sx = plotX + (rpmBin - xMin) / (xMax - xMin) * plotW;
float sy = plotY + (1 - (bin.torqueNm - yRightMin) / (yRightMax - yRightMin)) * plotH;
if (firstTorque) { torqueLine.Clear(); firstTorque = false; }
torqueLine.Append(new Vertex(new Vector2f(sx, sy), torqueColor));
}
else if (!firstTorque)
{
target.Draw(torqueLine);
torqueLine.Clear();
firstTorque = true;
}
}
if (!firstTorque) target.Draw(torqueLine);
if (currentRpm > 0 && currentRpm <= xMax && currentPowerKw > 0)
{
float sx = plotX + (currentRpm - xMin) / (xMax - xMin) * plotW;
float sy = plotY + (1 - (currentPowerKw - yLeftMin) / (yLeftMax - yLeftMin)) * plotH;
var dot = new CircleShape(2.5f)
{
FillColor = Color.White,
Position = new Vector2f(sx - 2.5f, sy - 2.5f)
};
target.Draw(dot);
}
}
// ---- Drawing helpers ----
protected Color PressureColor(float pressurePa)
{
float bar = pressurePa / 1e5f;
@@ -68,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;
@@ -107,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);
@@ -116,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;
}
@@ -157,5 +367,18 @@ namespace FluidSim.Tests
}
target.Draw(va);
}
protected void DrawLabel(RenderWindow target, string text, Vector2f position, Color fillColor, uint characterSize = 14)
{
if (Font == null) return;
var txt = new Text(Font)
{
DisplayedString = text,
Position = position,
FillColor = fillColor,
CharacterSize = characterSize
};
target.Draw(txt);
}
}
}

294
Scenarios/SingleCylScenario.cs
View File

@@ -1,6 +1,7 @@
using FluidSim.Components;
using FluidSim.Core;
using FluidSim.Interfaces;
using FluidSim.Utils;
using SFML.Graphics;
using SFML.System;
using System;
@@ -9,208 +10,200 @@ namespace FluidSim.Tests
{
public class SingleCylScenario : Scenario
{
// ---------- Engine components ----------
private Crankshaft crankshaft;
private Cylinder cylinder;
// ---------- Fluid network ----------
private PipeSystem pipeSystem;
private BoundarySystem boundaries;
private Solver solver;
// Volumes
private Volume0D intakePlenum;
// Ports
private Port plenumInlet, plenumOutlet;
private Volume0D exhaustCollector;
private Port colIn, colOut;
// Orifice / open‑end indices
private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
private int intakeOpenIdx, exhaustOpenIdx;
private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx;
private float[] orificeAreas;
private int intakeOpenIdx, exhaustOpenIdx;
// Sound
private SoundProcessor exhaustSound, intakeSound;
private OutdoorExhaustReverb reverb;
// ---------- Simulation state ----------
private double dt;
private int stepCount;
// ---------- Geometry (Lifan YX140) ----------
// Bore 56 mm, Stroke 57 mm, CR 9.5
private const float Bore = 0.056f;
private const float Stroke = 0.057f;
private const float ConRod = 0.110f; // typical for 57 mm stroke
private const float CompressionRatio = 9.5f;
// Valve diameters (intake 27 mm, exhaust 23 mm)
private const float IntakeValveDiam = 0.027f;
private const float ExhaustValveDiam = 0.023f;
private const float ValveLift = 0.006f; // 6 mm peak lift
// Valve timings (degrees, 720° four‑stroke)
// Intake: 15° BTDC → 45° ABDC
private const float IVO = 345f; // 15° BTDC
private const float IVC = 585f; // 45° ABDC (180°+45°)
// Exhaust: 45° BBDC → 15° ATDC
private const float EVO = 135f; // 45° BBDC (180°-45°)
private const float EVC = 375f; // 15° ATDC (360°+15°)
// Spark advance: 30° BTDC
private const float SparkAdv = 30f;
// Pipe / plenum sizes
private const float PipeDiam = 0.025f; // 25 mm intake / exhaust
private const float PipeArea = 0.00049087f; // π*D²/4
private const float PlenumVolume = 0.0005f; // 500 mL
private const float MaxThrottleArea = 1e-4f; // ~1 cm² (fully open)
// Pipe lengths and cell counts
private const float IntakeLenBefore = 0.15f; // 15 cm before throttle
private const float RunnerLen = 0.25f; // 25 cm runner
private const float ExhaustLen = 0.60f; // 60 cm exhaust
private const int CellsBefore = 6;
private const int CellsRunner = 10;
private const int CellsExhaust = 24;
private float _maxThrottleArea;
private float intakePipeArea, exhaustPipeArea;
private const float MaxBrakeTorque = 30.0f; // Nm at full load
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
// ---------- Crankshaft ----------
crankshaft = new Crankshaft(600); // start at ~600 RPM
crankshaft.Inertia = 0.2f;
crankshaft.FrictionConstant = 2.0f;
crankshaft.FrictionViscous = 0.04f;
// Throttle body diameter 44mm (typical for 250cc MX)
_maxThrottleArea = (float)Units.AreaFromDiameter(44 * Units.mm);
// ---------- Cylinder ----------
cylinder = new Cylinder(Bore, Stroke, ConRod, CompressionRatio,
IVO, IVC, EVO, EVC, crankshaft)
// ---- Crankshaft ----
crankshaft = new Crankshaft(2000);
crankshaft.Inertia = 0.02f; // kg·m² (crank + flywheel)
crankshaft.FrictionConstant = 3.0f; // Nm – bearings, rings, seals
crankshaft.FrictionViscous = 0.002f; // Nm/(rad/s) – oil windage
// ---- Cylinder (CRF250R) ----
float bore = 0.078f; // 78 mm
float stroke = 0.0522f; // 52.2 mm → 249.4 cc
float conRod = 0.1044f; // 2× stroke
float compRatio = 13.5f; // typical
// Valve events (high‑performance MX cam)
float ivo = 340f, ivc = 600f; // intake opens 20° BTDC (overlap), closes 60° ABDC
float evo = 120f, evc = 380f; // exhaust opens 60° BBDC, closes 20° ATDC
cylinder = new Cylinder(bore, stroke, conRod, compRatio,
ivo, ivc, evo, evc, crankshaft)
{
IntakeValveDiameter = IntakeValveDiam,
ExhaustValveDiameter = ExhaustValveDiam,
IntakeValveLift = ValveLift,
ExhaustValveLift = ValveLift,
SparkAdvance = SparkAdv,
EnergyVariationFraction = 0.03f, // small cycle‑to‑cycle variation
MisfireProbability = 0.0f
IntakeValveDiameter = 0.036f, // 36 mm
IntakeValveLift = 0.0095f, // 9.5 mm
ExhaustValveDiameter = 0.030f, // 30 mm
ExhaustValveLift = 0.0085f // 8.5 mm
};
// ---------- Pipe system ----------
int totalCells = CellsBefore + CellsRunner + CellsExhaust;
int[] pipeStart = { 0, CellsBefore, CellsBefore + CellsRunner };
int[] pipeEnd = { CellsBefore, CellsBefore + CellsRunner, totalCells };
// ---- Pipe system ----
int[] pipeStart = { 0, 10, 20 };
int[] pipeEnd = { 10, 20, 70 };
int totalCells = pipeEnd[^1];
float[] area = new float[totalCells];
float[] dx = new float[totalCells];
float[] areas = new float[totalCells];
float[] dxs = new float[totalCells];
float dxBefore = IntakeLenBefore / CellsBefore;
float dxRunner = RunnerLen / CellsRunner;
float dxExh = ExhaustLen / CellsExhaust;
float intakeDia = 0.040f; // 40 mm intake runner
float exhaustDia = 0.038f; // 38 mm exhaust primary
intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia;
exhaustPipeArea = MathF.PI * 0.25f * exhaustDia * exhaustDia;
float intakeLenBefore = 0.15f; // throttle body to plenum
float intakeLenRunner = 0.25f; // plenum to valve
float exhaustLen = 0.50f; // exhaust length
for (int i = 0; i < totalCells; i++)
{
areas[i] = PipeArea;
if (i < CellsBefore)
dxs[i] = dxBefore;
else if (i < CellsBefore + CellsRunner)
dxs[i] = dxRunner;
if (i < 10)
{
area[i] = intakePipeArea; dx[i] = intakeLenBefore / 10f;
}
else if (i < 20)
{
area[i] = intakePipeArea; dx[i] = intakeLenRunner / 10f;
}
else
dxs[i] = dxExh;
{
area[i] = exhaustPipeArea; dx[i] = exhaustLen / 50f;
}
}
float rho0 = 101325f / (287f * 300f);
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, areas, dxs,
rho0, 0f, 101325f);
pipeSystem.DampingMultiplier = 0.5f;
pipeSystem.EnergyRelaxationRate = 0f; // adiabatic pipes
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(PlenumVolume, 101325f, 300f);
plenumInlet = intakePlenum.CreatePort();
// ---- Volumes ----
intakePlenum = new Volume0D(1.0e-3f, 101325f, 300f); // 1 litre airbox
plenumInlet = intakePlenum.CreatePort();
plenumOutlet = intakePlenum.CreatePort();
exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // unused
colIn = exhaustCollector.CreatePort();
colOut = exhaustCollector.CreatePort();
// ---------- Boundary system ----------
// ---- Boundary system ----
boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
throttleAreaIdx = 0;
plenumRunnerIdx = 1;
intakeValveIdx = 2;
exhaustValveIdx = 3;
throttleAreaIdx = 0;
plenumRunnerAreaIdx = 1;
intakeValveIdx = 2;
exhaustValveIdx = 3;
// Open ends
boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, PipeArea);
// Open ends (pipe area = pipe cross‑section)
boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, intakePipeArea);
intakeOpenIdx = 0;
boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, PipeArea);
boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, exhaustPipeArea);
exhaustOpenIdx = 1;
// Orifices
// throttle – variable area, low discharge for restriction
boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false,
throttleAreaIdx, dischargeCoeff: 0.8f);
// plenum → runner
boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true,
plenumRunnerIdx, dischargeCoeff: 1.0f);
// intake valve
boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false,
intakeValveIdx, dischargeCoeff: 1.0f);
// exhaust valve
boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true,
exhaustValveIdx, dischargeCoeff: 1.0f);
boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.7f); // throttle
boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1.0f); // plenum→runner
boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1.0f); // intake valve
boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1.0f); // exhaust valve
orificeAreas = new float[4];
orificeAreas[plenumRunnerIdx] = PipeArea; // fixed full‑bore
orificeAreas[plenumRunnerAreaIdx] = intakePipeArea; // runner cross‑section (fixed)
// ---------- Solver ----------
solver = new Solver { SubStepCount = 5, EnableProfiling = false };
// ---- Solver ----
solver = new Solver { SubStepCount = 4, EnableProfiling = false };
solver.SetTimeStep(dt);
solver.SetPipeSystem(pipeSystem);
solver.SetBoundarySystem(boundaries);
solver.AddComponent(cylinder);
solver.AddComponent(intakePlenum);
solver.AddComponent(exhaustCollector);
// ---------- Sound ----------
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
// ---- Sound ----
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
reverb = new OutdoorExhaustReverb(sampleRate);
stepCount = 0;
Console.WriteLine("Single‑cylinder engine (YX140) ready.");
Console.WriteLine("CRF250R engine ready.");
}
public override float Process()
{
// ---- Crank and cylinder pre‑step ----
{
// Manual brake torque (0..30 Nm)
float loadTorque = Load * MaxBrakeTorque;
crankshaft.SetLoadTorque(loadTorque);
crankshaft.Step((float)dt);
cylinder.PreStep((float)dt);
// ---- Update variable areas ----
float throttledArea = MaxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1.0f);
float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
orificeAreas[throttleAreaIdx] = throttledArea;
orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
boundaries.SetOrificeAreas(orificeAreas);
// ---- Fluids step ----
solver.Step();
stepCount++;
// ---- Sound ----
float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
float exhaustDry = exhaustSound.Process(exhaustFlow);
float intakeDry = intakeSound.Process(intakeFlow);
if (stepCount % 2000 == 0)
if (stepCount % 1000 == 0)
{
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
Console.WriteLine($"Step {stepCount}, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar, " +
$"Throttle={Throttle * 100:F0}%");
float crankDeg = (crankshaft.CrankAngle + cylinder.PhaseOffset) * 180f / MathF.PI % 720f;
Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1}°, RPM={rpm:F0}, CylP={cylinder.Pressure/1e5f:F2} bar");
Console.WriteLine($" intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
var (r0L, u0L, p0L) = pipeSystem.GetInteriorStateLeft(0);
var (r0R, u0R, p0R) = pipeSystem.GetInteriorStateRight(0);
Console.WriteLine($" Pipe0 L: rho={r0L:F4} u={u0L:F3} p={p0L/1e5:F3}bar | R: rho={r0R:F4} u={u0R:F3} p={p0R/1e5:F3}bar");
var (r1L, u1L, p1L) = pipeSystem.GetInteriorStateLeft(1);
var (r1R, u1R, p1R) = pipeSystem.GetInteriorStateRight(1);
Console.WriteLine($" Pipe1 L: rho={r1L:F4} u={u1L:F3} p={p1L/1e5:F3}bar | R: rho={r1R:F4} u={u1R:F3} p={p1R/1e5:F3}bar");
var (r2L, u2L, p2L) = pipeSystem.GetInteriorStateLeft(2);
var (r2R, u2R, p2R) = pipeSystem.GetInteriorStateRight(2);
Console.WriteLine($" Pipe2 L: rho={r2L:F4} u={u2L:F3} p={p2L/1e5:F3}bar | R: rho={r2R:F4} u={u2R:F3} p={p2R/1e5:F3}bar");
Console.WriteLine($" Plenum P={intakePlenum.Pressure/1e5:F3}bar, mass={intakePlenum.Mass:E4} kg");
Console.WriteLine($" Cyl mass={cylinder.Mass:E4} kg");
}
return reverb.Process(exhaustDry + intakeDry);
return reverb.Process((intakeDry + exhaustDry) * 0.5f);
}
public override void Draw(RenderWindow target)
@@ -220,53 +213,54 @@ namespace FluidSim.Tests
float intakeY = winH / 2f - 40f;
float exhaustY = winH / 2f + 80f;
float leftX = 40f;
float openEndX = 40f;
// Intake open end marker
var om = new CircleShape(5f) { FillColor = Color.Cyan };
om.Position = new Vector2f(leftX - 5f, intakeY - 5f);
target.Draw(om);
float pipe1StartX = openEndX;
float pipe1EndX = pipe1StartX + 120f;
DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
// Pipe 0 – before throttle
float p0EndX = leftX + 80f;
DrawPipe(target, pipeSystem, 0, intakeY, leftX, p0EndX);
// Throttle symbol
float thrX = p0EndX + 5f;
var thr = new RectangleShape(new Vector2f(8f, 30f))
float throttleX = pipe1EndX + 5f;
var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
{
FillColor = Color.Yellow,
Position = new Vector2f(thrX, intakeY - 15f)
Position = new Vector2f(throttleX, intakeY - 15f)
};
target.Draw(thr);
target.Draw(throttleRect);
// Plenum volume
float plenW = 60f, plenH = 50f;
float plenLeftX = thrX + 12f;
float plenW = 60f, plenH = 80f;
float plenLeftX = throttleX + 10f;
float plenCenterX = plenLeftX + plenW / 2f;
float plenTopY = intakeY - plenH / 2f;
DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
// Pipe 1 – runner
float rStartX = plenLeftX + plenW + 10f;
float rEndX = rStartX + 100f;
DrawPipe(target, pipeSystem, 1, intakeY, rStartX, rEndX);
float runnerStartX = plenLeftX + plenW + 5f;
float runnerEndX = runnerStartX + 100f;
DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
// Cylinder
float cylCX = rEndX + 50f;
float cylCX = runnerEndX + 50f;
float cylTopY = intakeY - 120f;
float cylW = 80f, cylMaxH = 240f;
DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
// Pipe 2 – exhaust
float exhStartX = cylCX + cylW / 2f + 20f;
float exhEndX = winW - 60f;
DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
// Exhaust open end
var em = new CircleShape(5f) { FillColor = Color.Magenta };
em.Position = new Vector2f(exhEndX - 5f, exhaustY - 5f);
target.Draw(em);
// --- RPM & Power 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);
// --- Dyno curve ---
float torqueNm = crankshaft.AverageTorque;
UpdateDynoCurve(rpm, powerKw, torqueNm);
float graphX = winW - 410f;
float graphY = winH - 260f;
float graphW = 400f;
float graphH = 250f;
DrawDynoCurve(target, graphX, graphY, graphW, graphH, rpm, powerKw);
}
}
}

350
Scenarios/TwoStrokeScenario.cs Normal file
View File

@@ -0,0 +1,350 @@
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;
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 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;
cylinder = new TwoStrokeCylinder(bore, stroke, conRod, compRatio,
transferDuration, exhaustDuration,
crankshaft)
{
IntakeValveDiameter = 0.042f, // matched to intake pipe
IntakeValveLift = 0.015f,
ExhaustValveDiameter = 0.040f,
ExhaustValveLift = 0.013f
};
// ── 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
// --- 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
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];
// --- 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;
for (int i = 0; i < intakeCells; i++)
{ area[i] = intakePipeArea; dx[i] = intakeStubLen / intakeCells; }
for (int i = intakeCells; i < intakeCells + runnerCells; i++)
{ area[i] = intakePipeArea; dx[i] = intakeRunnerLen / runnerCells; }
// 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
exhaustHeaderArea = MathF.PI * 0.25f * headerDia * headerDia;
float bellyArea = MathF.PI * 0.25f * bellyDia * bellyDia;
float stingerArea = MathF.PI * 0.25f * stingerDia * stingerDia;
// 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 exhBase = intakeCells + runnerCells;
int idx = 0;
for (int i = exhBase; i < totalCells; i++, idx++)
{
if (idx < headerCells)
{
area[i] = exhaustHeaderArea;
dx[i] = headerLen / headerCells;
}
else if (idx < headerCells + diffuserCells)
{
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;
}
}
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
1.225f, 0f, 101325f);
pipeSystem.DampingMultiplier = 0.8f; // slightly less damping → stronger pulses
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);
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: 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);
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);
orificeAreas = new float[4];
orificeAreas[plenumRunnerIdx] = intakePipeArea; // runner always fully open
// ── Solver ────────────────────────────────────────────────────────────
// SubStepCount = 4 keeps CFL ≤ 1 for 5 mm cells at 44 100 Hz
solver = new Solver { SubStepCount = 4, EnableProfiling = false };
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 = 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}");
}
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);
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);
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");
}
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;
// Intake stub
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)
};
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 (expansion chamber)
float exhStartX = cylCX + 40f + 20f;
DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, winW - 60f, areaScale: 800f);
// HUD labels
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
float powerKw = crankshaft.AveragePower * 1e-3f;
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);
}
}
}