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

crank case based two stroke testing

Browse Source
This commit is contained in:
max
2026-06-10 00:13:14 +02:00
parent 56e9c2867a
commit 16498f8041
4 changed files with 406 additions and 308 deletions
Download Patch File Download Diff File Expand all files Collapse all files

150
Components/Crankcase.cs Normal file
View 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;
}
}
}
}
}

1
Components/Crankshaft.cs
View File

@@ -27,6 +27,7 @@ namespace FluidSim.Components
/// <summary>Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke.</summary>
public float CycleLength { get; set; } = 4f * MathF.PI;
public float CrankAngleRad => CrankAngle;
public Crankshaft(float initialRPM = 400f)
{

199
Components/TwoStrokeCylinder.cs
View File

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