Compare commits
4 Commits
aba9b76530
...
Testing
| Author | SHA1 | Date | |
|---|---|---|---|
| 16498f8041 | |||
| 56e9c2867a | |||
| 1240ebc33d | |||
| ac2eab6f83 |
150
Components/Crankcase.cs
Normal file
150
Components/Crankcase.cs
Normal file
@@ -0,0 +1,150 @@
|
||||
// ============================================================
|
||||
// File: Crankcase.cs
|
||||
// ============================================================
|
||||
using System.Collections.Generic;
|
||||
using FluidSim.Interfaces;
|
||||
|
||||
namespace FluidSim.Components
|
||||
{
|
||||
public class Crankcase : IComponent
|
||||
{
|
||||
private readonly Crankshaft _crankshaft;
|
||||
private readonly float _crankRadius, _conrodLength, _pistonArea;
|
||||
private readonly float _clearanceVolume, _obliquity;
|
||||
|
||||
private float _mass, _internalEnergy, _airFraction;
|
||||
public float Pressure { get; private set; }
|
||||
public float Temperature { get; private set; }
|
||||
public float Density => _mass / MathF.Max(Volume, 1e-12f);
|
||||
public float Volume { get; private set; }
|
||||
|
||||
public Port IntakePort { get; }
|
||||
public Port TransferPort { get; }
|
||||
private readonly List<Port> _ports;
|
||||
public IReadOnlyList<Port> Ports => _ports;
|
||||
|
||||
// FIX: store previous volume to calculate PdV work
|
||||
private float _prevVolume;
|
||||
|
||||
private const float Rgas = 287.0f;
|
||||
private const float Gamma = 1.4f;
|
||||
private const float Cv = Rgas / (Gamma - 1.0f);
|
||||
|
||||
public Crankcase(Crankshaft crankshaft,
|
||||
float crankRadius, float conrodLength, float bore,
|
||||
float clearanceVolume,
|
||||
float initialPressure, float initialTemperature)
|
||||
{
|
||||
_crankshaft = crankshaft;
|
||||
_crankRadius = crankRadius;
|
||||
_conrodLength = conrodLength;
|
||||
_pistonArea = MathF.PI * 0.25f * bore * bore;
|
||||
_clearanceVolume = clearanceVolume;
|
||||
_obliquity = crankRadius / conrodLength;
|
||||
|
||||
Pressure = initialPressure;
|
||||
Temperature = initialTemperature;
|
||||
float rho = initialPressure / (Rgas * initialTemperature);
|
||||
_mass = rho * clearanceVolume;
|
||||
_internalEnergy = _mass * Cv * initialTemperature;
|
||||
_airFraction = 1.0f;
|
||||
Volume = clearanceVolume;
|
||||
_prevVolume = Volume;
|
||||
|
||||
IntakePort = new Port { Owner = this };
|
||||
TransferPort = new Port { Owner = this };
|
||||
_ports = new List<Port> { IntakePort, TransferPort };
|
||||
}
|
||||
|
||||
public void PreStep(float dt)
|
||||
{
|
||||
// Save previous volume before updating
|
||||
_prevVolume = Volume;
|
||||
|
||||
float theta = _crankshaft.CrankAngleRad % (2f * MathF.PI);
|
||||
float cosTh = MathF.Cos(theta);
|
||||
float sinTh = MathF.Sin(theta);
|
||||
float term = MathF.Sqrt(1f - _obliquity * _obliquity * sinTh * sinTh);
|
||||
|
||||
// FIX: correct piston displacement: downstroke reduces crankcase volume
|
||||
float x = _crankRadius * (1f - cosTh) + _conrodLength * (1f - term);
|
||||
// Maximum volume at TDC (x = 0), minimum at BDC (x = stroke)
|
||||
float maxVolume = _clearanceVolume + _pistonArea * 2f * _crankRadius; // stroke = 2 * crankRadius
|
||||
Volume = maxVolume - _pistonArea * x;
|
||||
|
||||
// Update thermodynamic state using the new volume (before mass transfer)
|
||||
if (_mass > 1e-12f && Volume > 1e-12f)
|
||||
{
|
||||
Temperature = _internalEnergy / (_mass * Cv);
|
||||
Pressure = _mass * Rgas * Temperature / Volume;
|
||||
}
|
||||
}
|
||||
|
||||
public void UpdateState(float dt)
|
||||
{
|
||||
// ---- Mass and energy transport (identical to original) ----
|
||||
float mdotIn = IntakePort.MassFlowRate;
|
||||
float mdotOut = TransferPort.MassFlowRate;
|
||||
|
||||
float dm = (mdotIn - mdotOut) * dt;
|
||||
float dE = (mdotIn * IntakePort.SpecificEnthalpy
|
||||
- mdotOut * (Cv * Temperature + Pressure / MathF.Max(Density, 1e-12f))) * dt;
|
||||
float dY = (mdotIn * IntakePort.AirFraction - mdotOut * _airFraction) * dt;
|
||||
|
||||
_mass += dm;
|
||||
_internalEnergy += dE;
|
||||
if (_mass > 1e-12f)
|
||||
_airFraction = Math.Clamp((_airFraction * (_mass - dm) + dY) / _mass, 0f, 1f);
|
||||
|
||||
// ---- FIX: add mechanical work done BY the gas ON the piston ----
|
||||
// During a step, volume changed from _prevVolume to current Volume.
|
||||
// Work done BY gas = P * dV (if dV > 0, gas expands and does work, losing energy)
|
||||
float dV = Volume - _prevVolume;
|
||||
// Use average pressure during the step (approximate with current pressure)
|
||||
_internalEnergy -= Pressure * dV; // removes energy when volume increases
|
||||
|
||||
// Safety floors
|
||||
if (_mass < 1e-9f)
|
||||
{
|
||||
_mass = 1e-9f;
|
||||
_internalEnergy = _mass * Cv * 300f;
|
||||
_airFraction = 1f;
|
||||
}
|
||||
if (_internalEnergy < 0f)
|
||||
_internalEnergy = _mass * Cv * 300f;
|
||||
|
||||
// Final state update
|
||||
if (_mass > 1e-12f && Volume > 1e-12f)
|
||||
{
|
||||
Temperature = _internalEnergy / (_mass * Cv);
|
||||
Pressure = _mass * Rgas * Temperature / Volume;
|
||||
}
|
||||
else
|
||||
{
|
||||
Temperature = 300f;
|
||||
Pressure = 101325f;
|
||||
}
|
||||
|
||||
// Safety limits (unchanged, but now rarely triggered)
|
||||
const float safetyPressure = 1.0f;
|
||||
if (Pressure < safetyPressure && _mass > 1e-12f && Volume > 1e-12f)
|
||||
{
|
||||
Temperature = safetyPressure * Volume / (_mass * Rgas);
|
||||
_internalEnergy = _mass * Cv * Temperature;
|
||||
Pressure = safetyPressure;
|
||||
}
|
||||
|
||||
const float maxPressure = 5e5f;
|
||||
if (Pressure > maxPressure && _mass > 1e-12f && Volume > 1e-12f)
|
||||
{
|
||||
float targetMass = maxPressure * Volume / (Rgas * Temperature);
|
||||
if (_mass > targetMass)
|
||||
{
|
||||
_mass = targetMass;
|
||||
_internalEnergy = _mass * Cv * Temperature;
|
||||
Pressure = maxPressure;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -4,25 +4,52 @@ 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 float CrankAngleRad => CrankAngle;
|
||||
|
||||
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 +61,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;
|
||||
}
|
||||
|
||||
@@ -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.0f;
|
||||
public float CylinderWallArea = 0.02f;
|
||||
public float HeatTransferCoefficient = 100f;
|
||||
public float AmbientTemperature = 300f;
|
||||
|
||||
public float PhaseOffset; // rad
|
||||
|
||||
public float Volume => cylinderVolume;
|
||||
public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
|
||||
public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
|
||||
public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
|
||||
public float Mass => _airMass + _exhaustMass;
|
||||
public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
|
||||
public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
|
||||
|
||||
private float cylinderVolume, cylinderEnergy;
|
||||
private float _airMass, _exhaustMass;
|
||||
private float trappedAirMass, fuelMass, burnFraction;
|
||||
private bool combustionActive, fuelInjected;
|
||||
private float _energyFactor = 1f;
|
||||
private readonly Random _random = new Random();
|
||||
|
||||
private const float Gamma = 1.4f;
|
||||
private const float GasConstant = 287f;
|
||||
private const float MaxPressurePa = 200e5f;
|
||||
private const float MaxTemperatureK = 3500f;
|
||||
public override float IntakeValveArea =>
|
||||
MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
|
||||
public override float ExhaustValveArea =>
|
||||
MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
|
||||
|
||||
public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio,
|
||||
float ivo, float ivc, float evo, float evc, Crankshaft crankshaft)
|
||||
: base(bore, stroke, conRodLength, compressionRatio, crankshaft)
|
||||
{
|
||||
Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
|
||||
CompressionRatio = compressionRatio;
|
||||
IVO = ivo; IVC = ivc; EVO = evo; EVC = evc;
|
||||
Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
|
||||
|
||||
cylinderVolume = clearanceVolume;
|
||||
float initRho = 1.225f;
|
||||
_airMass = initRho * clearanceVolume;
|
||||
_exhaustMass = 0f;
|
||||
cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
|
||||
|
||||
IntakePort = new Port { Owner = this };
|
||||
ExhaustPort = new Port { Owner = this };
|
||||
_ports = new[] { IntakePort, ExhaustPort };
|
||||
}
|
||||
|
||||
private float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
|
||||
private float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
|
||||
private float CrankRadius => Stroke * 0.5f;
|
||||
private float Obliquity => CrankRadius / ConRodLength;
|
||||
|
||||
private float CrankDeg =>
|
||||
((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f;
|
||||
|
||||
public float ComputeVolume(float thetaRad)
|
||||
{
|
||||
float r = CrankRadius, l = ConRodLength;
|
||||
float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
|
||||
float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
|
||||
float x = r * (1f - cosTh) + l * (1f - term);
|
||||
float area = MathF.PI * 0.25f * Bore * Bore;
|
||||
return clearanceVolume + area * x;
|
||||
}
|
||||
|
||||
private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
|
||||
@@ -101,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
203
Components/EngineCylinder.cs
Normal 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
192
Components/TwoStrokeCylinder.cs
Normal file
192
Components/TwoStrokeCylinder.cs
Normal file
@@ -0,0 +1,192 @@
|
||||
// ============================================================
|
||||
// File: TwoStrokeCylinder.cs
|
||||
// ============================================================
|
||||
using System;
|
||||
using FluidSim.Interfaces;
|
||||
using FluidSim.Components; // for Crankcase (if in same namespace)
|
||||
|
||||
namespace FluidSim.Components
|
||||
{
|
||||
/// <summary>
|
||||
/// Two‑stroke cylinder with forced symmetrical port timings around BDC (180°).
|
||||
/// Uses crankcase back‑pressure for accurate pumping work.
|
||||
/// </summary>
|
||||
public class TwoStrokeCylinder : EngineCylinder
|
||||
{
|
||||
// --- Port timing (computed from durations) ---
|
||||
public float IVO => 180f - transferDuration / 2f;
|
||||
public float IVC => 180f + transferDuration / 2f;
|
||||
public float EVO => 180f - exhaustDuration / 2f;
|
||||
public float EVC => 180f + exhaustDuration / 2f;
|
||||
|
||||
private readonly float transferDuration; // degrees
|
||||
private readonly float exhaustDuration; // degrees
|
||||
|
||||
// --- Crankcase reference ---
|
||||
private Crankcase? _crankcase;
|
||||
|
||||
protected override float CycleLengthRad => 2f * MathF.PI;
|
||||
protected override float MaxCycleDeg => 360f;
|
||||
|
||||
public override float IntakeValveArea =>
|
||||
MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
|
||||
public override float ExhaustValveArea =>
|
||||
MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
|
||||
|
||||
public TwoStrokeCylinder(float bore, float stroke, float conRodLength,
|
||||
float compressionRatio,
|
||||
float transferDuration, float exhaustDuration,
|
||||
Crankshaft crankshaft)
|
||||
: base(bore, stroke, conRodLength, compressionRatio, crankshaft)
|
||||
{
|
||||
this.transferDuration = transferDuration;
|
||||
this.exhaustDuration = exhaustDuration;
|
||||
|
||||
if (EVO >= IVO)
|
||||
throw new ArgumentException("Exhaust must open before transfer port.");
|
||||
}
|
||||
|
||||
public void SetCrankcase(Crankcase crankcase)
|
||||
{
|
||||
_crankcase = crankcase;
|
||||
}
|
||||
|
||||
// ----- Valve lift -----
|
||||
private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
|
||||
{
|
||||
float deg = thetaDeg % 360f;
|
||||
if (deg < 0f) deg += 360f;
|
||||
|
||||
float effectiveOpen = opens;
|
||||
float effectiveClose = closes;
|
||||
if (closes < opens) effectiveClose += 360f;
|
||||
float duration = effectiveClose - effectiveOpen;
|
||||
if (duration <= 0f) return 0f;
|
||||
|
||||
float mapped = deg;
|
||||
if (mapped < opens) mapped += 360f;
|
||||
if (mapped < opens || mapped > effectiveClose) return 0f;
|
||||
|
||||
float rampDur = duration * 0.25f;
|
||||
float holdDur = duration - 2f * rampDur;
|
||||
|
||||
if (mapped >= opens && mapped < opens + rampDur)
|
||||
{
|
||||
float t = (mapped - opens) / rampDur;
|
||||
return peakLift * t * t * (3f - 2f * t);
|
||||
}
|
||||
else if (mapped >= opens + rampDur && mapped < opens + rampDur + holdDur)
|
||||
{
|
||||
return peakLift;
|
||||
}
|
||||
else if (mapped >= opens + rampDur + holdDur && mapped <= effectiveClose)
|
||||
{
|
||||
float t = (mapped - (opens + rampDur + holdDur)) / rampDur;
|
||||
return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
|
||||
}
|
||||
return 0f;
|
||||
}
|
||||
|
||||
protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
|
||||
{
|
||||
// Transfer port closing → fuel injection
|
||||
if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
|
||||
{
|
||||
trappedAirMass = _airMass;
|
||||
fuelMass = trappedAirMass / StoichiometricAFR;
|
||||
fuelInjected = true;
|
||||
}
|
||||
|
||||
// Spark every 360° at TDC (0°) minus advance
|
||||
float sparkAngle = (0f - SparkAdvance + 360f) % 360f;
|
||||
bool crossedSpark = false;
|
||||
if (prevDeg < sparkAngle && currDeg >= sparkAngle)
|
||||
crossedSpark = true;
|
||||
else if (prevDeg > sparkAngle && currDeg < sparkAngle)
|
||||
crossedSpark = true;
|
||||
|
||||
if (crossedSpark && !combustionActive && fuelInjected)
|
||||
{
|
||||
if (_random.NextDouble() < MisfireProbability)
|
||||
{
|
||||
combustionActive = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
combustionActive = true; burnFraction = 0f;
|
||||
float range = EnergyVariationFraction;
|
||||
_energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
|
||||
}
|
||||
}
|
||||
|
||||
if (combustionActive)
|
||||
{
|
||||
float angleSinceSpark = currDeg - sparkAngle;
|
||||
if (angleSinceSpark < 0f) angleSinceSpark += 360f;
|
||||
float newFraction = Wiebe(angleSinceSpark);
|
||||
if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance))
|
||||
{
|
||||
newFraction = 1f; combustionActive = false;
|
||||
float totalMass = _airMass + _exhaustMass;
|
||||
_airMass = 0f; _exhaustMass = totalMass;
|
||||
}
|
||||
fuelInjected = false;
|
||||
|
||||
float dFraction = newFraction - burnFraction;
|
||||
if (dFraction > 0f)
|
||||
{
|
||||
float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction;
|
||||
cylinderEnergy += dQ;
|
||||
_exhaustMass += fuelMass * dFraction;
|
||||
burnFraction = newFraction;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// ----- Override torque calculation to use crankcase back‑pressure -----
|
||||
public new void PreStep(float dt)
|
||||
{
|
||||
// Speed‑dependent spark advance
|
||||
float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
|
||||
SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);
|
||||
|
||||
float prevVolume = cylinderVolume;
|
||||
float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
|
||||
cylinderVolume = ComputeVolume(crankAngleRad);
|
||||
|
||||
float dV = cylinderVolume - prevVolume;
|
||||
|
||||
// Use crankcase pressure as back‑pressure, ambient if not set
|
||||
float backPressure = _crankcase?.Pressure ?? 101325f;
|
||||
float pRel = Pressure - backPressure;
|
||||
|
||||
float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
|
||||
float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
|
||||
float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
|
||||
float pistonArea = MathF.PI * 0.25f * Bore * Bore;
|
||||
Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
|
||||
|
||||
cylinderEnergy -= Pressure * dV;
|
||||
|
||||
float cycleLenDeg = 360f;
|
||||
float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % cycleLenDeg;
|
||||
float currDeg = crankAngleRad * 180f / MathF.PI % cycleLenDeg;
|
||||
|
||||
HandleCycleEvents(prevDeg, currDeg, dt);
|
||||
|
||||
// Heat loss
|
||||
float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
|
||||
(Temperature - AmbientTemperature) * dt;
|
||||
cylinderEnergy -= dQ_loss;
|
||||
|
||||
// Update port states
|
||||
float p = Pressure, rho = Density, T = Temperature;
|
||||
float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
|
||||
float af = AirFraction;
|
||||
IntakePort.Pressure = p; IntakePort.Density = rho;
|
||||
IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
|
||||
ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
|
||||
ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
|
||||
}
|
||||
}
|
||||
}
|
||||
166
Components/Vehicle.cs
Normal file
166
Components/Vehicle.cs
Normal 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;
|
||||
}
|
||||
}
|
||||
60
Program.cs
60
Program.cs
@@ -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;
|
||||
}
|
||||
}
|
||||
@@ -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);
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -32,53 +32,72 @@ namespace FluidSim.Tests
|
||||
private double dt;
|
||||
private int stepCount;
|
||||
|
||||
// Use a private field for the maximum throttle area, avoiding any base‑class conflicts
|
||||
private float _maxThrottleArea;
|
||||
|
||||
// pipe area for open end calculations
|
||||
private float pipeArea;
|
||||
private float intakePipeArea, exhaustPipeArea;
|
||||
private const float MaxBrakeTorque = 30.0f; // Nm at full load
|
||||
|
||||
public override void Initialize(int sampleRate)
|
||||
{
|
||||
dt = 1.0 / sampleRate;
|
||||
|
||||
// Maximum throttle area – independent of base class
|
||||
_maxThrottleArea = (float)Units.AreaFromDiameter(3 * Units.cm); // 1 cm²
|
||||
// Throttle body diameter 44mm (typical for 250cc MX)
|
||||
_maxThrottleArea = (float)Units.AreaFromDiameter(44 * Units.mm);
|
||||
|
||||
// ---- Crankshaft ----
|
||||
crankshaft = new Crankshaft(2000);
|
||||
crankshaft.Inertia = 0.01f;
|
||||
crankshaft.FrictionConstant = 2f;
|
||||
crankshaft.FrictionViscous = 0.0f;
|
||||
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 ----
|
||||
float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 11f;
|
||||
float ivo = 350f, ivc = 580f, evo = 120f, evc = 370f;
|
||||
cylinder = new Cylinder(bore, stroke, conRod, compRatio,
|
||||
ivo, ivc, evo, evc, crankshaft)
|
||||
{
|
||||
IntakeValveDiameter = 0.03f,
|
||||
IntakeValveLift = 0.005f,
|
||||
ExhaustValveDiameter = 0.028f,
|
||||
ExhaustValveLift = 0.005f
|
||||
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[] pipeStart = { 0, 10, 20 };
|
||||
int[] pipeEnd = { 10, 20, 70 };
|
||||
int totalCells = pipeEnd[^1]; // automatically 70, stays in sync
|
||||
int totalCells = pipeEnd[^1];
|
||||
float[] area = new float[totalCells];
|
||||
float[] dx = new float[totalCells];
|
||||
float pipeDiameter = 0.02f; // 2 cm
|
||||
pipeArea = MathF.PI * 0.25f * pipeDiameter * pipeDiameter;
|
||||
float areaVal = pipeArea;
|
||||
float intakeLenBefore = 0.2f, intakeLenRunner = 0.2f, exhaustLen = 0.4f;
|
||||
|
||||
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++)
|
||||
{
|
||||
area[i] = areaVal;
|
||||
if (i < 10) dx[i] = intakeLenBefore / 10f;
|
||||
else if (i < 20) dx[i] = intakeLenRunner / 10f;
|
||||
else dx[i] = exhaustLen / 50f;
|
||||
if (i < 10)
|
||||
{
|
||||
area[i] = intakePipeArea; dx[i] = intakeLenBefore / 10f;
|
||||
}
|
||||
else if (i < 20)
|
||||
{
|
||||
area[i] = intakePipeArea; dx[i] = intakeLenRunner / 10f;
|
||||
}
|
||||
else
|
||||
{
|
||||
area[i] = exhaustPipeArea; dx[i] = exhaustLen / 50f;
|
||||
}
|
||||
}
|
||||
|
||||
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
|
||||
@@ -88,10 +107,10 @@ namespace FluidSim.Tests
|
||||
pipeSystem.AmbientPressure = 101325f;
|
||||
|
||||
// ---- Volumes ----
|
||||
intakePlenum = new Volume0D(100e-6f, 101325f, 300f); // 100 mL
|
||||
intakePlenum = new Volume0D(1.0e-3f, 101325f, 300f); // 1 litre airbox
|
||||
plenumInlet = intakePlenum.CreatePort();
|
||||
plenumOutlet = intakePlenum.CreatePort();
|
||||
exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // 10 mL (unused but present)
|
||||
exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // unused
|
||||
colIn = exhaustCollector.CreatePort();
|
||||
colOut = exhaustCollector.CreatePort();
|
||||
|
||||
@@ -103,28 +122,20 @@ namespace FluidSim.Tests
|
||||
intakeValveIdx = 2;
|
||||
exhaustValveIdx = 3;
|
||||
|
||||
// Intake open end (pipe0 left)
|
||||
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;
|
||||
|
||||
// Throttle orifice (plenum inlet to pipe0 right)
|
||||
boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.2f);
|
||||
|
||||
// Plenum to runner (plenum outlet to pipe1 left)
|
||||
boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1f);
|
||||
|
||||
// Intake valve (cylinder intake to pipe1 right)
|
||||
boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1f);
|
||||
|
||||
// Exhaust valve (cylinder exhaust to pipe2 left)
|
||||
boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1f);
|
||||
|
||||
// Exhaust open end (pipe2 right)
|
||||
boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, pipeArea);
|
||||
boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, exhaustPipeArea);
|
||||
exhaustOpenIdx = 1;
|
||||
|
||||
// Orifices
|
||||
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[plenumRunnerAreaIdx] = areaVal; // fixed plenum->runner area
|
||||
orificeAreas[plenumRunnerAreaIdx] = intakePipeArea; // runner cross‑section (fixed)
|
||||
|
||||
// ---- Solver ----
|
||||
solver = new Solver { SubStepCount = 4, EnableProfiling = false };
|
||||
@@ -136,22 +147,26 @@ namespace FluidSim.Tests
|
||||
solver.AddComponent(exhaustCollector);
|
||||
|
||||
// ---- Sound ----
|
||||
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
|
||||
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
|
||||
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
|
||||
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
|
||||
reverb = new OutdoorExhaustReverb(sampleRate);
|
||||
|
||||
stepCount = 0;
|
||||
Console.WriteLine("TestScenario ready.");
|
||||
Console.WriteLine("CRF250R engine ready.");
|
||||
}
|
||||
|
||||
public override float Process()
|
||||
{
|
||||
// Manual brake torque (0..30 Nm)
|
||||
float loadTorque = Load * MaxBrakeTorque;
|
||||
crankshaft.SetLoadTorque(loadTorque);
|
||||
|
||||
crankshaft.Step((float)dt);
|
||||
cylinder.PreStep((float)dt);
|
||||
|
||||
// Update variable orifice areas – use the private _maxThrottleArea
|
||||
float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1f);
|
||||
float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
|
||||
orificeAreas[throttleAreaIdx] = throttledArea;
|
||||
|
||||
orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
|
||||
orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
|
||||
boundaries.SetOrificeAreas(orificeAreas);
|
||||
@@ -159,7 +174,6 @@ namespace FluidSim.Tests
|
||||
solver.Step();
|
||||
stepCount++;
|
||||
|
||||
// Retrieve open‑end mass flows for sound synthesis
|
||||
float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
|
||||
float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
|
||||
|
||||
@@ -169,31 +183,27 @@ namespace FluidSim.Tests
|
||||
if (stepCount % 1000 == 0)
|
||||
{
|
||||
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
|
||||
float crankDeg = crankshaft.CrankAngle; // degrees (0–720)
|
||||
Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1} deg, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar");
|
||||
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}");
|
||||
|
||||
// Pipe 0 (intake before throttle)
|
||||
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");
|
||||
|
||||
// Pipe 1 (runner)
|
||||
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");
|
||||
|
||||
// Pipe 2 (exhaust)
|
||||
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");
|
||||
|
||||
// Plenum and cylinder mass
|
||||
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(intakeDry + exhaustDry);
|
||||
return reverb.Process((intakeDry + exhaustDry) * 0.5f);
|
||||
}
|
||||
|
||||
public override void Draw(RenderWindow target)
|
||||
@@ -205,12 +215,10 @@ namespace FluidSim.Tests
|
||||
float exhaustY = winH / 2f + 80f;
|
||||
float openEndX = 40f;
|
||||
|
||||
// Intake pipe before throttle (pipe 0)
|
||||
float pipe1StartX = openEndX;
|
||||
float pipe1EndX = pipe1StartX + 120f;
|
||||
DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
|
||||
|
||||
// Throttle symbol
|
||||
float throttleX = pipe1EndX + 5f;
|
||||
var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
|
||||
{
|
||||
@@ -219,28 +227,40 @@ namespace FluidSim.Tests
|
||||
};
|
||||
target.Draw(throttleRect);
|
||||
|
||||
// Plenum
|
||||
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);
|
||||
|
||||
// Runner pipe (pipe 1)
|
||||
float runnerStartX = plenLeftX + plenW + 5f;
|
||||
float runnerEndX = runnerStartX + 100f;
|
||||
DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
|
||||
|
||||
// Cylinder
|
||||
float cylCX = runnerEndX + 50f;
|
||||
float cylTopY = intakeY - 120f;
|
||||
float cylW = 80f, cylMaxH = 240f;
|
||||
DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
|
||||
|
||||
// Exhaust pipe (pipe 2)
|
||||
float exhStartX = cylCX + cylW / 2f + 20f;
|
||||
float exhEndX = winW - 60f;
|
||||
DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
|
||||
|
||||
// --- 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
288
Scenarios/TwoStrokeScenario.cs
Normal file
288
Scenarios/TwoStrokeScenario.cs
Normal file
@@ -0,0 +1,288 @@
|
||||
using FluidSim.Components;
|
||||
using FluidSim.Core;
|
||||
using FluidSim.Interfaces;
|
||||
using FluidSim.Utils;
|
||||
using SFML.Graphics;
|
||||
using SFML.System;
|
||||
using System;
|
||||
|
||||
namespace FluidSim.Tests
|
||||
{
|
||||
public class TwoStrokeScenario : Scenario
|
||||
{
|
||||
private Crankshaft crankshaft;
|
||||
private TwoStrokeCylinder cylinder;
|
||||
private Crankcase crankcase;
|
||||
|
||||
private PipeSystem pipeSystem;
|
||||
private BoundarySystem boundaries;
|
||||
private Solver solver;
|
||||
|
||||
private Volume0D intakePlenum;
|
||||
private Port plenumInlet, plenumOutlet;
|
||||
|
||||
private int throttleAreaIdx, reedInletIdx, reedOutletIdx,
|
||||
transferInletIdx, transferOutletIdx, exhaustValveIdx;
|
||||
private float[] orificeAreas;
|
||||
private int intakeOpenIdx, exhaustOpenIdx;
|
||||
|
||||
private SoundProcessor exhaustSound, intakeSound;
|
||||
private OutdoorExhaustReverb reverb;
|
||||
|
||||
private double dt;
|
||||
private int stepCount;
|
||||
|
||||
private float maxThrottleArea;
|
||||
private float intakePipeArea, reedPipeArea, transferPipeArea, exhaustHeaderArea;
|
||||
private bool reedOpen;
|
||||
|
||||
public override void Initialize(int sampleRate)
|
||||
{
|
||||
dt = 1.0 / sampleRate;
|
||||
|
||||
maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm);
|
||||
|
||||
// ---- Crankshaft ----
|
||||
crankshaft = new Crankshaft(3000);
|
||||
crankshaft.CycleLength = 2f * MathF.PI;
|
||||
crankshaft.Inertia = 0.01f;
|
||||
crankshaft.FrictionConstant = 1.0f;
|
||||
crankshaft.FrictionViscous = 0.002f;
|
||||
|
||||
// ---- Cylinder (125cc) ----
|
||||
float bore = 0.054f, stroke = 0.0545f, conRod = 0.110f, compRatio = 7.2f;
|
||||
float transferDur = 140f, exhaustDur = 195f;
|
||||
cylinder = new TwoStrokeCylinder(bore, stroke, conRod, compRatio,
|
||||
transferDur, exhaustDur, crankshaft)
|
||||
{
|
||||
// FIX: realistic transfer port diameter (was 40mm)
|
||||
IntakeValveDiameter = 0.030f, // 30 mm
|
||||
IntakeValveLift = 0.010f,
|
||||
ExhaustValveDiameter = 0.040f,
|
||||
ExhaustValveLift = 0.010f
|
||||
};
|
||||
|
||||
// ---- Crankcase ----
|
||||
float crankRadius = stroke * 0.5f;
|
||||
float ccClearance = 150e-6f;
|
||||
crankcase = new Crankcase(crankshaft, crankRadius, conRod, bore,
|
||||
ccClearance, 101325f, 300f);
|
||||
cylinder.SetCrankcase(crankcase);
|
||||
|
||||
// ---- Pipe system ----
|
||||
int intakeCells = 8;
|
||||
int reedCells = 4;
|
||||
int transferCells = 8;
|
||||
int exhaustCells = 60;
|
||||
int totalCells = intakeCells + reedCells + transferCells + exhaustCells;
|
||||
|
||||
int[] pipeStart = {
|
||||
0,
|
||||
intakeCells,
|
||||
intakeCells + reedCells,
|
||||
intakeCells + reedCells + transferCells
|
||||
};
|
||||
int[] pipeEnd = {
|
||||
intakeCells,
|
||||
intakeCells + reedCells,
|
||||
intakeCells + reedCells + transferCells,
|
||||
totalCells
|
||||
};
|
||||
|
||||
float[] area = new float[totalCells];
|
||||
float[] dx = new float[totalCells];
|
||||
|
||||
float intakeDia = 0.042f, reedDia = 0.040f, transferDia = 0.040f;
|
||||
intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia;
|
||||
reedPipeArea = MathF.PI * 0.25f * reedDia * reedDia;
|
||||
transferPipeArea = MathF.PI * 0.25f * transferDia * transferDia;
|
||||
|
||||
for (int i = 0; i < intakeCells; i++)
|
||||
{ area[i] = intakePipeArea; dx[i] = 0.100f / intakeCells; }
|
||||
|
||||
for (int i = intakeCells; i < intakeCells + reedCells; i++)
|
||||
{ area[i] = reedPipeArea; dx[i] = 0.030f / reedCells; }
|
||||
|
||||
for (int i = intakeCells + reedCells; i < intakeCells + reedCells + transferCells; i++)
|
||||
{ area[i] = transferPipeArea; dx[i] = 0.200f / transferCells; }
|
||||
|
||||
float hdrD = 0.040f, hdrL = 0.130f;
|
||||
float difEndD = 0.070f, difL = 0.250f;
|
||||
float belL = 0.220f;
|
||||
float convEndD = 0.028f, convL = 0.160f;
|
||||
float stiL = 0.080f;
|
||||
float totL = hdrL + difL + belL + convL + stiL;
|
||||
exhaustHeaderArea = MathF.PI * 0.25f * hdrD * hdrD;
|
||||
float bellyArea = MathF.PI * 0.25f * difEndD * difEndD;
|
||||
float stingerArea = MathF.PI * 0.25f * convEndD * convEndD;
|
||||
|
||||
int exhStart = intakeCells + reedCells + transferCells;
|
||||
int hdrC = (int)(exhaustCells * hdrL / totL);
|
||||
int difC = (int)(exhaustCells * difL / totL);
|
||||
int belC = (int)(exhaustCells * belL / totL);
|
||||
int conC = (int)(exhaustCells * convL / totL);
|
||||
int stiC = exhaustCells - hdrC - difC - belC - conC;
|
||||
|
||||
int idx = 0;
|
||||
for (int i = exhStart; i < totalCells; i++)
|
||||
{
|
||||
if (idx < hdrC) { area[i] = exhaustHeaderArea; dx[i] = hdrL / hdrC; }
|
||||
else if (idx < hdrC + difC)
|
||||
{
|
||||
float t = (idx - hdrC) / (float)(difC - 1);
|
||||
float dia = hdrD + (difEndD - hdrD) * t;
|
||||
area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = difL / difC;
|
||||
}
|
||||
else if (idx < hdrC + difC + belC) { area[i] = bellyArea; dx[i] = belL / belC; }
|
||||
else if (idx < hdrC + difC + belC + conC)
|
||||
{
|
||||
float t = (idx - hdrC - difC - belC) / (float)(conC - 1);
|
||||
float dia = difEndD + (convEndD - difEndD) * t;
|
||||
area[i] = MathF.PI * 0.25f * dia * dia; dx[i] = convL / conC;
|
||||
}
|
||||
else { area[i] = stingerArea; dx[i] = stiL / stiC; }
|
||||
idx++;
|
||||
}
|
||||
|
||||
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
|
||||
1.225f, 0f, 101325f);
|
||||
pipeSystem.DampingMultiplier = 0.8f;
|
||||
pipeSystem.EnergyRelaxationRate = 0.4f;
|
||||
|
||||
// ---- Volumes ----
|
||||
intakePlenum = new Volume0D(0.5e-3f, 101325f, 300f);
|
||||
plenumInlet = intakePlenum.CreatePort();
|
||||
plenumOutlet = intakePlenum.CreatePort();
|
||||
|
||||
// ---- Boundary system ----
|
||||
boundaries = new BoundarySystem(pipeSystem, maxOrifices: 6, maxOpenEnds: 2);
|
||||
throttleAreaIdx = 0;
|
||||
reedInletIdx = 1;
|
||||
reedOutletIdx = 2;
|
||||
transferInletIdx = 3;
|
||||
transferOutletIdx = 4;
|
||||
exhaustValveIdx = 5;
|
||||
|
||||
boundaries.AddOpenEnd(0, true, 101325f, intakePipeArea);
|
||||
intakeOpenIdx = 0;
|
||||
boundaries.AddOpenEnd(3, false, 101325f, stingerArea);
|
||||
exhaustOpenIdx = 1;
|
||||
|
||||
boundaries.AddOrifice(plenumInlet, 0, false, throttleAreaIdx, 0.72f);
|
||||
boundaries.AddOrifice(plenumOutlet, 1, true, reedInletIdx, 1.0f);
|
||||
boundaries.AddOrifice(crankcase.IntakePort, 1, false, reedOutletIdx, 0.9f);
|
||||
boundaries.AddOrifice(crankcase.TransferPort,2, true, transferInletIdx,1.0f);
|
||||
boundaries.AddOrifice(cylinder.IntakePort, 2, false, transferOutletIdx,1.0f);
|
||||
boundaries.AddOrifice(cylinder.ExhaustPort, 3, true, exhaustValveIdx, 0.7f);
|
||||
|
||||
orificeAreas = new float[6];
|
||||
orificeAreas[reedInletIdx] = reedPipeArea;
|
||||
orificeAreas[reedOutletIdx] = 0f;
|
||||
orificeAreas[transferInletIdx] = transferPipeArea;
|
||||
orificeAreas[transferOutletIdx] = 0f;
|
||||
|
||||
// ---- Solver ----
|
||||
solver = new Solver { SubStepCount = 4 };
|
||||
solver.SetTimeStep(dt);
|
||||
solver.SetPipeSystem(pipeSystem);
|
||||
solver.SetBoundarySystem(boundaries);
|
||||
solver.AddComponent(cylinder);
|
||||
solver.AddComponent(crankcase);
|
||||
solver.AddComponent(intakePlenum);
|
||||
|
||||
// ---- Sound ----
|
||||
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
|
||||
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
|
||||
reverb = new OutdoorExhaustReverb(sampleRate);
|
||||
|
||||
stepCount = 0;
|
||||
Console.WriteLine("Two‑Stroke engine ready.");
|
||||
}
|
||||
|
||||
public override float Process()
|
||||
{
|
||||
const float reedMargin = 200f;
|
||||
if (crankcase.Pressure < intakePlenum.Pressure - reedMargin)
|
||||
reedOpen = true;
|
||||
else if (crankcase.Pressure > intakePlenum.Pressure + reedMargin)
|
||||
reedOpen = false;
|
||||
|
||||
float throttledFraction = Throttle;
|
||||
if (throttledFraction < 0.001f) throttledFraction = 0f;
|
||||
throttledFraction = Math.Clamp(throttledFraction, 0f, 1f);
|
||||
float throttledArea = maxThrottleArea * throttledFraction;
|
||||
|
||||
orificeAreas[throttleAreaIdx] = throttledArea;
|
||||
orificeAreas[reedOutletIdx] = reedOpen ? reedPipeArea : 0f;
|
||||
orificeAreas[transferOutletIdx] = cylinder.IntakeValveArea;
|
||||
orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
|
||||
boundaries.SetOrificeAreas(orificeAreas);
|
||||
|
||||
if (stepCount < 20000)
|
||||
crankshaft.AddTorque(5.0f);
|
||||
|
||||
// FIX: update crankshaft BEFORE volumes, so crankcase and cylinder see the same new angle
|
||||
crankshaft.Step((float)dt);
|
||||
cylinder.PreStep((float)dt);
|
||||
crankcase.PreStep((float)dt);
|
||||
|
||||
solver.Step();
|
||||
stepCount++;
|
||||
|
||||
float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
|
||||
float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
|
||||
|
||||
float exhaustDry = exhaustSound.Process(exhaustFlow);
|
||||
float intakeDry = intakeSound.Process(intakeFlow);
|
||||
|
||||
if (stepCount % 2000 == 0)
|
||||
{
|
||||
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
|
||||
Console.WriteLine($"Step {stepCount} | RPM={rpm:F0} | CylP={cylinder.Pressure/1e5f:F2} bar | CCP={crankcase.Pressure/1e5f:F3} bar | Plenum={intakePlenum.Pressure/1e5f:F3} bar | Reed={reedOpen}");
|
||||
}
|
||||
|
||||
return reverb.Process((intakeDry + exhaustDry) * 0.5f);
|
||||
}
|
||||
|
||||
public override void Draw(RenderWindow target)
|
||||
{
|
||||
float winW = target.GetView().Size.X;
|
||||
float winH = target.GetView().Size.Y;
|
||||
|
||||
float startX = 40f;
|
||||
float endX = winW - 80f;
|
||||
|
||||
DrawPipe(target, pipeSystem, 0, winH * 0.25f, startX, startX + 120f);
|
||||
var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
|
||||
{
|
||||
FillColor = Color.Yellow,
|
||||
Position = new Vector2f(startX + 125f, winH * 0.25f - 15f)
|
||||
};
|
||||
target.Draw(throttleRect);
|
||||
float plenX = startX + 140f;
|
||||
DrawVolume(target, intakePlenum, plenX + 30f, winH * 0.25f - 25f, 60f, 50f);
|
||||
|
||||
float reedStartX = plenX + 70f;
|
||||
DrawPipe(target, pipeSystem, 1, winH * 0.25f, reedStartX, reedStartX + 30f);
|
||||
|
||||
float transStartX = reedStartX + 40f;
|
||||
DrawPipe(target, pipeSystem, 2, winH * 0.45f, transStartX, transStartX + 120f);
|
||||
|
||||
float cylCX = transStartX + 180f;
|
||||
float cylTopY = winH * 0.45f - 90f;
|
||||
DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f);
|
||||
|
||||
float exhStartX = cylCX + 60f;
|
||||
DrawPipe(target, pipeSystem, 3, winH * 0.65f, exhStartX, endX, areaScale: 800f);
|
||||
|
||||
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
|
||||
float powerKw = crankshaft.AveragePower * 1e-3f;
|
||||
DrawLabel(target, $"RPM: {rpm:F0}", new Vector2f(20, 90), Color.White, 24);
|
||||
DrawLabel(target, $"Power: {powerKw:F2} kW", new Vector2f(20, 115), Color.White, 24);
|
||||
|
||||
float torqueNm = crankshaft.AverageTorque;
|
||||
UpdateDynoCurve(rpm, powerKw, torqueNm);
|
||||
DrawDynoCurve(target, winW - 410f, winH - 260f, 400f, 250f, rpm, powerKw);
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user