diff --git a/Components/Crankshaft.cs b/Components/Crankshaft.cs
index 99ea43e..240bb4b 100644
--- a/Components/Crankshaft.cs
+++ b/Components/Crankshaft.cs
@@ -7,7 +7,7 @@ namespace FluidSim.Components
{
public double AngularVelocity { get; set; } // rad/s
public double CrankAngle { get; set; } // rad, 0 … 4π (four‑stroke cycle)
- public double PreviousAngle { get; private set; } // for TDC detection
+ public double PreviousAngle { get; set; } // ← now has public setter
public double Inertia { get; set; } = 0.2;
public double FrictionConstant { get; set; } = 2.0; // N·m
@@ -15,7 +15,6 @@ namespace FluidSim.Components
private double externalTorque;
- /// Idle speed before any combustion torque is applied.
public Crankshaft(double initialRPM = 400.0)
{
AngularVelocity = initialRPM * 2.0 * Math.PI / 60.0;
@@ -27,20 +26,23 @@ namespace FluidSim.Components
public void Step(double dt)
{
- // Save previous angle
+ // Catch NaN before it propagates
+ if (double.IsNaN(AngularVelocity) || double.IsInfinity(AngularVelocity))
+ AngularVelocity = 0.0;
+ if (double.IsNaN(externalTorque) || double.IsInfinity(externalTorque))
+ externalTorque = 0.0;
+
PreviousAngle = CrankAngle;
- // Friction
double friction = FrictionConstant * Math.Sign(AngularVelocity) + FrictionViscous * AngularVelocity;
double netTorque = externalTorque - friction;
double alpha = netTorque / Inertia;
AngularVelocity += alpha * dt;
- if (AngularVelocity < 0) AngularVelocity = 0; // stall
+ if (AngularVelocity < 0) AngularVelocity = 0;
CrankAngle += AngularVelocity * dt;
- // Wrap to [0, 4π)
if (CrankAngle >= 4.0 * Math.PI)
CrankAngle -= 4.0 * Math.PI;
else if (CrankAngle < 0)
diff --git a/Components/EngineCylinder.cs b/Components/EngineCylinder.cs
index aa322fe..645fa29 100644
--- a/Components/EngineCylinder.cs
+++ b/Components/EngineCylinder.cs
@@ -1,4 +1,3 @@
-// EngineCylinder.cs (in Core namespace)
using System;
using FluidSim.Components;
@@ -11,21 +10,32 @@ namespace FluidSim.Core
private double bore, stroke, conRodLength, compressionRatio;
private double pistonArea;
- private double V_disp, V_clear;
- private double maxOrificeArea;
- private double valveOpenStart = 120.0 * Math.PI / 180.0;
- private double valveOpenEnd = 480.0 * Math.PI / 180.0;
- private double valveRampWidth = 30.0 * Math.PI / 180.0;
- public double OrificeArea => ValveLift() * maxOrificeArea;
+ public double V_disp { get; private set; }
+ public double V_clear { get; private set; }
+ public bool ignition = false;
+ // ---- Exhaust valve ----
+ private double exhMaxOrificeArea;
+ private double exhValveOpenStart = 120.0 * Math.PI / 180.0; // 120° (EVO)
+ private double exhValveOpenEnd = 480.0 * Math.PI / 180.0; // 480° (EVC)
+ private double exhValveRampWidth = 30.0 * Math.PI / 180.0;
+ public double ExhaustOrificeArea => ExhaustValveLift() * exhMaxOrificeArea;
+
+ // ---- Intake valve ----
+ private double intMaxOrificeArea;
+ private double intValveOpenStart = 380.0 * Math.PI / 180.0; // 380° (IVO)
+ private double intValveOpenEnd = 560.0 * Math.PI / 180.0; // 560° (IVC)
+ private double intValveRampWidth = 30.0 * Math.PI / 180.0;
+ public double IntakeOrificeArea => IntakeValveLift() * intMaxOrificeArea;
+
+ // ---- Combustion ----
public double TargetPeakPressure { get; set; } = 50.0 * 101325.0;
private const double PeakTemperature = 2500.0;
private bool burnInProgress = false;
- private double burnStartAngle; // full cycle angle when ignition began
+ private double burnStartAngle; // cycle angle (0–4π)
private double burnDuration = 40.0 * Math.PI / 180.0;
private double targetBurnEnergy;
- private double totalBurnMass;
private double preIgnitionMass, preIgnitionInternalEnergy;
private Random rand = new Random();
@@ -35,46 +45,64 @@ namespace FluidSim.Core
public int CombustionCount { get; private set; }
public int MisfireCount { get; private set; }
+ // Cycle‑aware angle (0 – 4π)
+ public double CycleAngle => crankshaft.CrankAngle % (4.0 * Math.PI);
+ private double prevCycleAngle;
+
+ // Piston position fraction (0 = TDC, 1 = BDC)
+ public double PistonPositionFraction
+ {
+ get
+ {
+ double currentVol = Cylinder.Volume;
+ if (currentVol <= V_clear) return 0.0;
+ if (currentVol >= V_clear + V_disp) return 1.0;
+ return (currentVol - V_clear) / V_disp;
+ }
+ }
+
public EngineCylinder(Crankshaft crankshaft,
double bore, double stroke, double compressionRatio,
- double pipeArea, int sampleRate)
+ double exhPipeArea, double intPipeArea, int sampleRate)
{
this.crankshaft = crankshaft;
this.bore = bore;
this.stroke = stroke;
conRodLength = 2.0 * stroke;
this.compressionRatio = compressionRatio;
- maxOrificeArea = pipeArea;
+ exhMaxOrificeArea = exhPipeArea;
+ intMaxOrificeArea = intPipeArea;
pistonArea = Math.PI / 4.0 * bore * bore;
V_disp = pistonArea * stroke;
V_clear = V_disp / (compressionRatio - 1.0);
- // Initial compressed charge at TDC (no burn)
- double T_bdc = 300.0;
- double p_bdc = 101325.0;
+ // Start at BDC with fresh ambient charge
double V_bdc = V_clear + V_disp;
- double freshMass = p_bdc * V_bdc / (287.0 * T_bdc);
- double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0);
- double p_tdc = p_bdc * Math.Pow(V_bdc / V_clear, 1.4);
+ double p_amb = 101325.0;
+ double T_amb = 300.0;
+ double rho0 = p_amb / (287.0 * T_amb);
+ double mass0 = rho0 * V_bdc;
+ double energy0 = p_amb * V_bdc / (1.4 - 1.0);
- Cylinder = new Volume0D(V_clear, p_tdc, T_bdc * Math.Pow(V_bdc / V_clear, 1.4 - 1.0), sampleRate)
+ Cylinder = new Volume0D(V_bdc, p_amb, T_amb, sampleRate)
{
Gamma = 1.4,
GasConstant = 287.0
};
- Cylinder.Volume = V_clear;
- Cylinder.Mass = freshMass;
- Cylinder.InternalEnergy = p_tdc * V_clear / (1.4 - 1.0);
+ Cylinder.Volume = V_bdc;
+ Cylinder.Mass = mass0;
+ Cylinder.InternalEnergy = energy0;
+
+ prevCycleAngle = CycleAngle;
preIgnitionMass = Cylinder.Mass;
preIgnitionInternalEnergy = Cylinder.InternalEnergy;
}
- // ---- Piston kinematics (uses full cycle angle for position) ----
+ // ---- Piston kinematics ----
private (double volume, double dvdt) PistonKinematics(double cycleAngle)
{
- // Slider-crank uses 0–2π, but we want the same motion for 0–2π (power/exhaust) and 2π–4π (intake/compression)
double theta = cycleAngle % (2.0 * Math.PI);
double R = stroke / 2.0;
double cosT = Math.Cos(theta);
@@ -90,26 +118,34 @@ namespace FluidSim.Core
return (V, dvdt);
}
- // ---- Valve lift ----
- private double ValveLift()
+ // ---- Valve lifts (cycle‑aware) ----
+ private double ExhaustValveLift()
{
- double cycleRad = crankshaft.CrankAngle;
- if (cycleRad < valveOpenStart || cycleRad > valveOpenEnd)
- return 0.0;
-
- double duration = valveOpenEnd - valveOpenStart;
- double ramp = valveRampWidth;
- double t = (cycleRad - valveOpenStart) / duration;
+ double a = CycleAngle;
+ if (a < exhValveOpenStart || a > exhValveOpenEnd) return 0.0;
+ double duration = exhValveOpenEnd - exhValveOpenStart;
+ double ramp = exhValveRampWidth;
+ double t = (a - exhValveOpenStart) / duration;
double rampFrac = ramp / duration;
-
- if (t < rampFrac)
- return t / rampFrac;
- else if (t > 1.0 - rampFrac)
- return (1.0 - t) / rampFrac;
- else
- return 1.0;
+ if (t < rampFrac) return t / rampFrac;
+ if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
+ return 1.0;
}
+ private double IntakeValveLift()
+ {
+ double a = CycleAngle;
+ if (a < intValveOpenStart || a > intValveOpenEnd) return 0.0;
+ double duration = intValveOpenEnd - intValveOpenStart;
+ double ramp = intValveRampWidth;
+ double t = (a - intValveOpenStart) / duration;
+ double rampFrac = ramp / duration;
+ if (t < rampFrac) return t / rampFrac;
+ if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
+ return 1.0;
+ }
+
+ // ---- Wiebe burn fraction ----
private double WiebeFraction(double angleFromIgnition)
{
if (angleFromIgnition >= burnDuration) return 1.0;
@@ -137,33 +173,24 @@ namespace FluidSim.Core
return force * lever;
}
+ // ---- TDC detection (power stroke, at angle 0 mod 4π) ----
+ private bool DetectTDCPowerStroke()
+ {
+ double current = CycleAngle;
+ double previous = prevCycleAngle;
+ prevCycleAngle = current;
+ return (previous > 3.8 * Math.PI && current < 0.2 * Math.PI);
+ }
+
public void Step(double dt)
{
- double cycleAngle = crankshaft.CrankAngle;
- double prevAngle = crankshaft.PreviousAngle;
-
- // ----- TDC crossing detection (power stroke) -----
- // Power stroke TDC occurs at angle 0 (mod 4π). We detect when PreviousAngle was near 4π and CrankAngle wraps to near 0.
- bool crossingTDC = (prevAngle > 3.8 * Math.PI && cycleAngle < 0.2 * Math.PI) // normal forward
- || (prevAngle < 0.2 * Math.PI && cycleAngle > 3.8 * Math.PI); // (rare backward, ignore)
+ bool crossingTDC = DetectTDCPowerStroke();
if (crossingTDC)
{
misfireCurrent = rand.NextDouble() < MisfireProbability;
- // Fresh charge: trapped at BDC, compressed isentropically to V_clear
- double T_bdc = 300.0;
- double p_bdc = 101325.0;
- double V_bdc = V_clear + V_disp;
- double freshMass = p_bdc * V_bdc / (287.0 * T_bdc);
- double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0);
- double gamma = 1.4;
- double p_tdc = p_bdc * Math.Pow(V_bdc / V_clear, gamma);
-
- Cylinder.Volume = V_clear;
- Cylinder.Mass = freshMass;
- Cylinder.InternalEnergy = p_tdc * V_clear / (gamma - 1.0);
-
+ // *** Always capture the state at TDC, whether we burn or not ***
preIgnitionMass = Cylinder.Mass;
preIgnitionInternalEnergy = Cylinder.InternalEnergy;
@@ -171,43 +198,48 @@ namespace FluidSim.Core
{
MisfireCount++;
}
- else
+ else if (ignition)
{
- double V = V_clear;
- targetBurnEnergy = TargetPeakPressure * V / (gamma - 1.0);
- totalBurnMass = TargetPeakPressure * V / (287.0 * PeakTemperature);
+ double V = Cylinder.Volume;
+ targetBurnEnergy = TargetPeakPressure * V / (Cylinder.Gamma - 1.0);
+ if (double.IsNaN(targetBurnEnergy))
+ targetBurnEnergy = 101325.0 * V / (Cylinder.Gamma - 1.0);
burnInProgress = true;
- burnStartAngle = cycleAngle;
+ burnStartAngle = CycleAngle;
CombustionCount++;
}
}
- // ----- Burn progress -----
if (burnInProgress)
{
- double angleFromIgnition = cycleAngle - burnStartAngle;
- if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI; // wrap if needed
+ double angleFromIgnition = CycleAngle - burnStartAngle;
+ if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI;
if (angleFromIgnition >= burnDuration)
{
- Cylinder.Mass = totalBurnMass;
Cylinder.InternalEnergy = targetBurnEnergy;
burnInProgress = false;
}
else
{
double fraction = WiebeFraction(angleFromIgnition);
- Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction) + targetBurnEnergy * fraction;
- Cylinder.Mass = preIgnitionMass * (1.0 - fraction) + totalBurnMass * fraction;
+ Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction)
+ + targetBurnEnergy * fraction;
+ Cylinder.Mass = preIgnitionMass;
}
}
- // ----- Piston motion -----
- var (vol, dvdt) = PistonKinematics(cycleAngle);
+ var (vol, dvdt) = PistonKinematics(CycleAngle);
Cylinder.Volume = vol;
Cylinder.Dvdt = dvdt;
- // ----- Torque contribution -----
+ if (double.IsNaN(Cylinder.Pressure) || double.IsNaN(Cylinder.Temperature) || Cylinder.Mass < 1e-9)
+ {
+ double V = Math.Max(vol, V_clear);
+ Cylinder.Mass = 1.225 * V;
+ Cylinder.InternalEnergy = 101325.0 * V / (1.4 - 1.0);
+ }
+
double torque = ComputeTorque();
crankshaft.AddTorque(torque);
}
diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
index 6c414cb..8f20679 100644
--- a/Components/Pipe1D.cs
+++ b/Components/Pipe1D.cs
@@ -106,6 +106,12 @@ namespace FluidSim.Components
public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p;
+ public float GetFaceMassFlux(int faceIndex)
+ {
+ if (faceIndex < 0 || faceIndex > _n) return 0f;
+ return _fluxM[faceIndex];
+ }
+
public void SetGhostLeft(double rho, double u, double p)
{
_rhoGhostL = (float)rho;
diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs
index 6fee6d9..a5b83ff 100644
--- a/Components/Volume0D.cs
+++ b/Components/Volume0D.cs
@@ -1,4 +1,4 @@
-using System;
+using System;
namespace FluidSim.Components
{
@@ -15,10 +15,10 @@ namespace FluidSim.Components
private double _dt;
- public double Density => Mass / Volume;
- public double Pressure => (Gamma - 1.0) * InternalEnergy / Volume;
- public double Temperature => Pressure / (Density * GasConstant);
- public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
+ public double Density => Mass / Math.Max(Volume, 1e-12);
+ public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12);
+ public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
+ public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12);
public double MassFlowRateIn { get; set; }
public double SpecificEnthalpyIn { get; set; }
@@ -46,18 +46,20 @@ namespace FluidSim.Components
Mass += dm;
InternalEnergy += dE;
- // Safety: if mass becomes extremely small, reset internal energy to zero
- if (Mass < 1e-12)
+ // ---- ABSOLUTE SAFEGUARD: keep at least 1 µg of gas at ambient pressure ----
+ double minMass = 1e-9;
+ double V = Math.Max(Volume, 1e-12);
+ if (Mass < minMass)
{
- Mass = 0.0;
- InternalEnergy = 0.0;
+ Mass = minMass;
+ InternalEnergy = 5000.0 * V / (Gamma - 1.0); // 0.05 bar, not ambient
}
- else if (InternalEnergy < 1e-12)
+ else if (InternalEnergy < 0.0)
{
- InternalEnergy = 0.0;
+ InternalEnergy = 101325.0 * V / (Gamma - 1.0);
}
- // Avoid negative mass/energy
+ // Final non‑negative clamp
if (Mass < 0.0) Mass = 0.0;
if (InternalEnergy < 0.0) InternalEnergy = 0.0;
}
diff --git a/Core/Constants.cs b/Core/Constants.cs
new file mode 100644
index 0000000..ca5d277
--- /dev/null
+++ b/Core/Constants.cs
@@ -0,0 +1,11 @@
+namespace FluidSim.Core
+{
+ public static class Constants
+ {
+ public const double Gamma = 1.4;
+ public const double R_gas = 287.0; // J/(kg·K)
+ public const double P_amb = 101325.0; // Pa
+ public const double T_amb = 300.0; // K
+ public static readonly double Rho_amb = P_amb / (R_gas * T_amb); // ≈ 1.177 kg/m³
+ }
+}
\ No newline at end of file
diff --git a/Core/NozzleFlow.cs b/Core/NozzleFlow.cs
index 5b40733..904cf00 100644
--- a/Core/NozzleFlow.cs
+++ b/Core/NozzleFlow.cs
@@ -9,7 +9,6 @@ namespace FluidSim.Core
out double massFlow, out double rhoFace, out double uFace, out double pFace,
double gamma = 1.4)
{
- // Default fallback (no flow)
massFlow = 0.0;
rhoFace = 0.0;
uFace = 0.0;
@@ -29,16 +28,43 @@ namespace FluidSim.Core
double pr = downstreamPressure / p0;
double choked = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
- if (pr < choked) pr = choked;
- double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
- if (double.IsNaN(M)) return;
+ // If pr > 1, flow is INTO the cylinder (reverse), so we swap the roles.
+ bool reverse = (pr > 1.0);
+ if (reverse)
+ {
+ // Treat the cylinder as the downstream, the pipe as the upstream.
+ double p_up = downstreamPressure;
+ double T_up = 300.0; // pipe temperature (ambient)
+ double rho_up = downstreamPressure / (R * T_up);
- uFace = M * Math.Sqrt(gamma * R * T0);
- rhoFace = rho0 * Math.Pow(pr, 1.0 / gamma);
- pFace = p0 * pr;
+ double pr_rev = p0 / p_up; // now cylinder / pipe
+ if (pr_rev < choked) pr_rev = choked;
+
+ double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr_rev, -(gamma - 1.0) / gamma) - 1.0));
+ if (double.IsNaN(M)) return;
+
+ // Flow from pipe INTO cylinder (positive mass flow into volume)
+ uFace = M * Math.Sqrt(gamma * R * T_up);
+ rhoFace = rho_up * Math.Pow(pr_rev, 1.0 / gamma);
+ pFace = p_up * pr_rev;
+ massFlow = rhoFace * uFace * area;
+ // massFlow is positive = into cylinder
+ }
+ else
+ {
+ // Normal flow out of cylinder
+ if (pr < choked) pr = choked;
+
+ double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
+ if (double.IsNaN(M)) return;
+
+ uFace = M * Math.Sqrt(gamma * R * T0);
+ rhoFace = rho0 * Math.Pow(pr, 1.0 / gamma);
+ pFace = p0 * pr;
+ massFlow = -rhoFace * uFace * area; // negative = out of cylinder
+ }
- massFlow = rhoFace * uFace * area;
if (double.IsNaN(massFlow) || double.IsInfinity(massFlow))
massFlow = 0.0;
}
diff --git a/Core/OutdoorExhaustReverb.cs b/Core/OutdoorExhaustReverb.cs
index 3712771..23dc997 100644
--- a/Core/OutdoorExhaustReverb.cs
+++ b/Core/OutdoorExhaustReverb.cs
@@ -4,95 +4,124 @@ namespace FluidSim.Core
{
public class OutdoorExhaustReverb
{
- // ---- Geometry ----
- private const float GroundReflDelay = 0.008f; // 8 ms (≈1.3 m)
- private const float WallRefl1Delay = 0.045f; // ≈15 m
- private const float WallRefl2Delay = 0.080f; // ≈27 m
+ // ========== Early reflection delays (stereo: left/right) ==========
+ private readonly DelayLine groundL, groundR;
+ private readonly DelayLine wall1L, wall1R;
+ private readonly DelayLine wall2L, wall2R;
- private DelayLine groundRefl;
- private DelayLine wallRefl1;
- private DelayLine wallRefl2;
-
- // ---- FDN for late diffuse tail ----
- private const int FDN_CHANNELS = 8; // dense, realistic
- private DelayLine[] fdnDelays;
- private float[] fdnState;
- private OrthonormalMixer mixer; // energy‑preserving mixing
- private LowPassFilter[] channelFilters; // per‑channel air absorption
+ // ========== Diffuse tail FDNs (left/right each with 8 channels) ==========
+ private const int FDN_CHANNELS = 8;
+ private readonly DelayLine[] fdnL, fdnR;
+ private readonly float[] stateL, stateR;
+ private readonly OrthonormalMixer mixerL, mixerR;
+ private readonly LowPassFilter[] filterL, filterR;
public float DryMix { get; set; } = 1.0f;
public float EarlyMix { get; set; } = 0.5f;
public float TailMix { get; set; } = 0.9f;
- public float Feedback { get; set; } = 0.75f; // safe range 0.7‑0.9
- public float DampingFreq { get; set; } = 6000f; // Hz, above which air absorbs strongly
- public float MatrixCoeff { get; set; } = 0.5f; // (kept for compatibility, not used)
+ public float Feedback { get; set; } = 0.75f; // safe range 0.7‑0.9
+ public float DampingFreq { get; set; } = 6000f; // Hz
public OutdoorExhaustReverb(int sampleRate)
{
- // Early reflection lines
- groundRefl = new DelayLine((int)(sampleRate * GroundReflDelay));
- wallRefl1 = new DelayLine((int)(sampleRate * WallRefl1Delay));
- wallRefl2 = new DelayLine((int)(sampleRate * WallRefl2Delay));
+ // Early reflections – left/right offset by ~1‑2 ms for stereo width
+ groundL = new DelayLine((int)(sampleRate * 0.008)); // 8 ms
+ groundR = new DelayLine((int)(sampleRate * 0.010)); // 10 ms
+ wall1L = new DelayLine((int)(sampleRate * 0.045));
+ wall1R = new DelayLine((int)(sampleRate * 0.047));
+ wall2L = new DelayLine((int)(sampleRate * 0.080));
+ wall2R = new DelayLine((int)(sampleRate * 0.082));
- // FDN delays: prime numbers for dense modal density (70‑150 ms)
- int[] baseLengths = { 3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009 };
- fdnDelays = new DelayLine[FDN_CHANNELS];
+ // FDN delay lengths – prime numbers, offset between L/R
+ int[] lengthsL = { 3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009 };
+ int[] lengthsR = { 3217, 4049, 5531, 7043, 8537, 10037, 11519, 13033 };
+ fdnL = new DelayLine[FDN_CHANNELS];
+ fdnR = new DelayLine[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
{
- int len = Math.Min(baseLengths[i], (int)(sampleRate * 0.25)); // max 250 ms
- fdnDelays[i] = new DelayLine(len);
+ int lenL = Math.Min(lengthsL[i], (int)(sampleRate * 0.25));
+ int lenR = Math.Min(lengthsR[i], (int)(sampleRate * 0.25));
+ fdnL[i] = new DelayLine(lenL);
+ fdnR[i] = new DelayLine(lenR);
}
- fdnState = new float[FDN_CHANNELS];
- mixer = new OrthonormalMixer(FDN_CHANNELS);
+ stateL = new float[FDN_CHANNELS];
+ stateR = new float[FDN_CHANNELS];
+ mixerL = new OrthonormalMixer(FDN_CHANNELS);
+ mixerR = new OrthonormalMixer(FDN_CHANNELS);
- // Air absorption: a gentle first‑order low‑pass per channel
- channelFilters = new LowPassFilter[FDN_CHANNELS];
- float initialCutoff = DampingFreq;
+ filterL = new LowPassFilter[FDN_CHANNELS];
+ filterR = new LowPassFilter[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
- channelFilters[i] = new LowPassFilter(sampleRate, initialCutoff);
+ {
+ filterL[i] = new LowPassFilter(sampleRate, DampingFreq);
+ filterR[i] = new LowPassFilter(sampleRate, DampingFreq);
+ }
}
- public float Process(float drySample)
+ /// Stereo reverb – returns (left, right) sample pair.
+ public (float left, float right) ProcessStereo(float drySample)
{
// ---- Early reflections ----
- float g = groundRefl.ReadWrite(drySample * 0.8f);
- float w1 = wallRefl1.ReadWrite(drySample * 0.5f);
- float w2 = wallRefl2.ReadWrite(drySample * 0.4f);
- float early = (g + w1 + w2) * EarlyMix;
+ float gL = groundL.ReadWrite(drySample * 0.8f);
+ float gR = groundR.ReadWrite(drySample * 0.8f);
+ float w1L = wall1L.ReadWrite(drySample * 0.5f);
+ float w1R = wall1R.ReadWrite(drySample * 0.5f);
+ float w2L = wall2L.ReadWrite(drySample * 0.4f);
+ float w2R = wall2R.ReadWrite(drySample * 0.4f);
- // ---- FDN diffuse tail ----
- // Read the delayed outputs (which were stored last iteration)
- float[] delOut = new float[FDN_CHANNELS];
+ float earlyL = (gL + w1L + w2L) * EarlyMix;
+ float earlyR = (gR + w1R + w2R) * EarlyMix;
+
+ // ---- Read diffuse tail ----
+ float[] delOutL = new float[FDN_CHANNELS];
+ float[] delOutR = new float[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
- delOut[i] = fdnDelays[i].Read();
+ {
+ delOutL[i] = fdnL[i].Read();
+ delOutR[i] = fdnR[i].Read();
+ }
- // Mix the delayed outputs with the orthonormal matrix -> scattered signals
- mixer.Process(delOut, fdnState); // result written into fdnState
+ // Mix via orthonormal matrix
+ float[] mixL = new float[FDN_CHANNELS];
+ float[] mixR = new float[FDN_CHANNELS];
+ mixerL.Process(delOutL, mixL);
+ mixerR.Process(delOutR, mixR);
- // Add fresh input to all channels
+ // Feedback + air absorption
for (int i = 0; i < FDN_CHANNELS; i++)
- fdnState[i] = drySample * 0.15f + Feedback * fdnState[i];
+ {
+ stateL[i] = drySample * 0.15f + Feedback * mixL[i];
+ stateL[i] = filterL[i].Process(stateL[i]);
+ fdnL[i].Write(stateL[i]);
- // Air absorption: per‑channel one‑pole low‑pass
+ stateR[i] = drySample * 0.15f + Feedback * mixR[i];
+ stateR[i] = filterR[i].Process(stateR[i]);
+ fdnR[i].Write(stateR[i]);
+ }
+
+ float tailL = 0.0f, tailR = 0.0f;
for (int i = 0; i < FDN_CHANNELS; i++)
- fdnState[i] = channelFilters[i].Process(fdnState[i]);
+ {
+ tailL += delOutL[i];
+ tailR += delOutR[i];
+ }
+ tailL *= TailMix;
+ tailR *= TailMix;
- // Write the new states into the delay lines
- for (int i = 0; i < FDN_CHANNELS; i++)
- fdnDelays[i].Write(fdnState[i]);
-
- // The tail output is the sum of the delayed outputs *before* the loop
- float tailSum = 0f;
- for (int i = 0; i < FDN_CHANNELS; i++)
- tailSum += delOut[i];
- float tail = tailSum * TailMix;
-
- // Final mix
- return drySample * DryMix + early + tail;
+ float left = drySample * DryMix + earlyL + tailL;
+ float right = drySample * DryMix + earlyR + tailR;
+ return (left, right);
}
- // ---------- Helper classes (same as before but with separate Read/Write) ----------
+ /// Mono fallback – sums left+right / 2.
+ public float Process(float drySample)
+ {
+ var (l, r) = ProcessStereo(drySample);
+ return (l + r) * 0.5f;
+ }
+
+ // ========== Helper classes ==========
private class DelayLine
{
private float[] buffer;
@@ -100,19 +129,13 @@ namespace FluidSim.Core
public DelayLine(int length)
{
buffer = new float[Math.Max(length, 1)];
- writePos = 0;
- }
- // Separated Read/Write to avoid ringing with immediate feedback
- public float Read()
- {
- return buffer[writePos];
}
+ public float Read() => buffer[writePos];
public void Write(float value)
{
buffer[writePos] = value;
writePos = (writePos + 1) % buffer.Length;
}
- // Old combined method (not used in FDN, only for early reflections)
public float ReadWrite(float value)
{
float outVal = buffer[writePos];
@@ -124,8 +147,7 @@ namespace FluidSim.Core
private class LowPassFilter
{
- private float b0, a1;
- private float y1;
+ private float b0, a1, y1;
private float sampleRate;
public LowPassFilter(int sampleRate, float cutoff)
{
@@ -137,7 +159,7 @@ namespace FluidSim.Core
float w = 2 * (float)Math.PI * cutoff / sampleRate;
float a0 = 1 + w;
b0 = w / a0;
- a1 = (1 - w) / a0; // first‑order low‑pass
+ a1 = (1 - w) / a0;
}
public float Process(float x)
{
@@ -147,18 +169,13 @@ namespace FluidSim.Core
}
}
- ///
- /// Computes a fast orthonormal mixing matrix (like Hadamard, but energy‑preserving).
- ///
private class OrthonormalMixer
{
private int size;
- public OrthonormalMixer(int size) { this.size = size; }
+ public OrthonormalMixer(int size) => this.size = size;
public void Process(float[] input, float[] output)
{
- // Simple energy‑conserving “allpass” mixing:
- // Use a Householder reflection: y = (2/n) * sum(x) * ones - x
float sum = 0;
for (int i = 0; i < size; i++) sum += input[i];
float factor = 2.0f / size;
diff --git a/Core/PipeVolumeConnection.cs b/Core/PipeVolumeConnection.cs
index 40eb6f7..4ff4b74 100644
--- a/Core/PipeVolumeConnection.cs
+++ b/Core/PipeVolumeConnection.cs
@@ -9,6 +9,8 @@ namespace FluidSim.Core
public bool IsPipeLeftEnd { get; }
public double OrificeArea { get; set; }
+ public double LastMassFlowIntoVolume { get; set; }
+
public PipeVolumeConnection(Volume0D vol, Pipe1D pipe, bool isPipeLeftEnd, double orificeArea)
{
Volume = vol;
diff --git a/Core/Solver.cs b/Core/Solver.cs
index e2efaad..66c019d 100644
--- a/Core/Solver.cs
+++ b/Core/Solver.cs
@@ -35,52 +35,85 @@ namespace FluidSim.Core
public float Step()
{
- // 1. Compute nozzle flows and update volumes (once per audio sample)
+ // 1. For each connection, handle flow or closed wall
foreach (var conn in _connections)
{
+ double area = conn.OrificeArea;
+ if (area < 1e-12) // valve closed → treat as solid wall
+ {
+ conn.Volume.MassFlowRateIn = 0.0;
+ conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy; // not used
+
+ // Set ghost to a reflective wall (u = -u_pipe, same p, ρ)
+ int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1;
+ double rho = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6);
+ double p = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0);
+ double u = conn.Pipe.GetCellVelocity(cellIdx);
+ if (conn.IsPipeLeftEnd)
+ conn.Pipe.SetGhostLeft(rho, -u, p);
+ else
+ conn.Pipe.SetGhostRight(rho, -u, p);
+ continue;
+ }
+
+ // Valve open → use the nozzle model
double downstreamPressure = conn.IsPipeLeftEnd
? conn.Pipe.GetCellPressure(0)
: conn.Pipe.GetCellPressure(conn.Pipe.GetCellCount() - 1);
- NozzleFlow.Compute(conn.Volume, conn.OrificeArea, downstreamPressure,
+ NozzleFlow.Compute(conn.Volume, area, downstreamPressure,
out double mdot, out double rhoFace, out double uFace, out double pFace,
gamma: conn.Volume.Gamma);
- // Limit mass flow to available mass
+ // Clamp mdot to available mass
double maxMdot = conn.Volume.Mass / _dt;
+ conn.LastMassFlowIntoVolume = mdot;
if (mdot > maxMdot) mdot = maxMdot;
if (mdot < -maxMdot) mdot = -maxMdot;
- conn.Volume.MassFlowRateIn = -mdot;
- conn.Volume.SpecificEnthalpyIn = (conn.Volume.Gamma / (conn.Volume.Gamma - 1.0)) *
- (conn.Volume.Pressure / Math.Max(conn.Volume.Density, 1e-12));
+ conn.Volume.MassFlowRateIn = mdot;
+ // enthalpy: if inflow, use pipe enthalpy; if outflow, use cylinder enthalpy
+ if (mdot >= 0)
+ {
+ int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1;
+ double pPipe = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0);
+ double rhoPipe = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6);
+ conn.Volume.SpecificEnthalpyIn = (conn.Volume.Gamma / (conn.Volume.Gamma - 1.0)) * pPipe / rhoPipe;
+ }
+ else
+ {
+ conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy;
+ }
+
+ // Integrate the volume
conn.Volume.Integrate(_dt);
+ // Set ghost from nozzle face state (but don't allow zero density)
+ if (rhoFace < 1e-6) rhoFace = Constants.Rho_amb;
+ if (pFace < 100.0) pFace = Constants.P_amb;
if (conn.IsPipeLeftEnd)
conn.Pipe.SetGhostLeft(rhoFace, uFace, pFace);
else
conn.Pipe.SetGhostRight(rhoFace, uFace, pFace);
}
- // 2. Determine required sub‑steps
+ // 2. Sub‑step pipes
int nSub = 1;
foreach (var p in _pipes)
nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
double dtSub = _dt / nSub;
- // 3. Sub‑step loop for pipes
for (int sub = 0; sub < nSub; sub++)
foreach (var p in _pipes)
p.SimulateSingleStep(dtSub);
- // 4. Clear ghost flags
+ // 3. Clear ghost flags
foreach (var p in _pipes)
p.ClearGhostFlag();
- // 5. Return raw mass flow from the first pipe’s open end (assumed exhaust tailpipe)
+ // 4. Return exhaust tailpipe mass flow
if (_pipes.Count > 0)
return (float)_pipes[0].GetOpenEndMassFlow();
-
return 0f;
}
}
diff --git a/Program.cs b/Program.cs
index 201c0cc..278e544 100644
--- a/Program.cs
+++ b/Program.cs
@@ -25,7 +25,7 @@ public class Program
// Throttle smoothing
private static double targetThrottle = 0.0; // 1.0 when W is pressed, 0.0 otherwise
private static double currentThrottle = 0.0;
- private const double ThrottleSmoothing = 8.0; // rate of change
+ private const double ThrottleSmoothing = 40.0; // rate of change
private static volatile bool running = true;
@@ -39,7 +39,7 @@ public class Program
window.KeyPressed += OnKeyPressed;
var soundEngine = new SoundEngine(bufferCapacity: 16384);
- soundEngine.Volume = 70;
+ soundEngine.Volume = 100;
soundEngine.Start();
scenario = new EngineScenario();
diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs
index ceda9d2..c0f513f 100644
--- a/Scenarios/EngineScenario.cs
+++ b/Scenarios/EngineScenario.cs
@@ -13,139 +13,197 @@ namespace FluidSim.Core
private Crankshaft crankshaft;
private EngineCylinder engineCyl;
private Pipe1D exhaustPipe;
- private PipeVolumeConnection coupling;
- private SoundProcessor soundProcessor;
+ private Pipe1D intakePipe;
+ private PipeVolumeConnection couplingExhaust;
+ private PipeVolumeConnection couplingIntake;
+ private SoundProcessor exhaustSoundProcessor;
+ private SoundProcessor intakeSoundProcessor;
private OutdoorExhaustReverb reverb;
- private Port exitPort = new Port();
+ private Port exhaustPort = new Port();
+ private Port intakePort = new Port();
private double dt;
- private double pipeArea;
+ private double exhPipeArea, intPipeArea;
private const double AmbientPressure = 101325.0;
private double time;
private int stepCount = 0;
- private const int LogInterval = 10000;
+ private const int LogInterval = 1000;
- // Throttle 0..1
- public double Throttle { get; set; } = 0.0; // start with a light idle throttle
-
- // ---- Realistic combustion parameters ----
- private const double FullLoadPeakPressure = 70.0 * 101325.0; // 15 bar
-
- // ---- Idle speed governor ----
- private const double TargetIdleRPM = 800.0; // rad/s = RPM * π/30, we'll convert
+ public double Throttle { get; set; } = 0.15;
+ private const double FullLoadPeakPressure = 60.0 * Units.bar;
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
- // ---- Crankshaft: inertia + friction that gives ~800 RPM at idle ----
- crankshaft = new Crankshaft(initialRPM: 600.0) // start a bit low
+ // Crankshaft
+ crankshaft = new Crankshaft(initialRPM: 2000.0)
{
- Inertia = 0.005, // slightly heavier flywheel
- FrictionConstant = 0.8, // static friction
- FrictionViscous = 0.01 // viscous (linear with RPM)
+ Inertia = 0.05,
+ FrictionConstant = 0.5,
+ FrictionViscous = 0.01
};
- // ---- Pipe: add a tiny bit of damping to prevent unrealistic shocks ----
- double pipeLength = 2;
- double pipeRadius = 0.1;
- pipeArea = Math.PI * pipeRadius * pipeRadius;
- exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 60);
+ // Exhaust pipe (longer, larger)
+ double exhLength = 1;
+ double exhRadius = 0.02;
+ exhPipeArea = Math.PI * exhRadius * exhRadius;
+ exhaustPipe = new Pipe1D(exhLength, exhPipeArea, sampleRate, forcedCellCount: 100);
exhaustPipe.SetUniformState(1.225, 0.0, AmbientPressure);
- exhaustPipe.DampingMultiplier = 5;
- exhaustPipe.EnergyRelaxationRate = 50;
+ exhaustPipe.DampingMultiplier = 0.0;
+ exhaustPipe.EnergyRelaxationRate = 100.0f;
- // ---- Cylinder ----
+ // Intake pipe (shorter, narrower)
+ double intLength = 1;
+ double intRadius = 0.01;
+ intPipeArea = Math.PI * intRadius * intRadius;
+ intakePipe = new Pipe1D(intLength, intPipeArea, sampleRate, forcedCellCount: 50);
+ intakePipe.SetUniformState(1.225, 0.0, AmbientPressure);
+
+ // Cylinder (starts at BDC, fresh charge)
engineCyl = new EngineCylinder(crankshaft,
- bore: 0.065, stroke: 0.0565, compressionRatio: 10.0,
- pipeArea: pipeArea, sampleRate: sampleRate);
+ bore: 0.065, stroke: 0.0565, compressionRatio: 8.0,
+ exhPipeArea: exhPipeArea, intPipeArea: intPipeArea, sampleRate: sampleRate);
+ engineCyl.ignition = true;
- // ---- Coupling ----
- coupling = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0);
+ // Set crank to BDC (180°) and sync
+ crankshaft.CrankAngle = Math.PI;
+ crankshaft.PreviousAngle = Math.PI; // make sure this property is settable (public setter)
- // ---- Solver ----
+ // Couplings
+ couplingExhaust = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0);
+ couplingIntake = new PipeVolumeConnection(engineCyl.Cylinder, intakePipe, false, orificeArea: 0.0);
+
+ // Solver
solver = new Solver();
solver.SetTimeStep(dt);
solver.AddVolume(engineCyl.Cylinder);
solver.AddPipe(exhaustPipe);
- solver.AddConnection(coupling);
+ solver.AddPipe(intakePipe);
+ solver.AddConnection(couplingExhaust);
+ solver.AddConnection(couplingIntake);
solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, AmbientPressure);
+ solver.SetPipeBoundary(intakePipe, true, BoundaryType.GhostCell); // cylinder side – left
+ solver.SetPipeBoundary(intakePipe, false, BoundaryType.OpenEnd, AmbientPressure); // ambient side – right
- // ---- Sound processor (stable version) ----
- soundProcessor = new SoundProcessor(sampleRate, pipeRadius * 2);
- soundProcessor.Gain = 0.00001f;
+ // Sound
+ exhaustSoundProcessor = new SoundProcessor(sampleRate, exhRadius * 2);
+ exhaustSoundProcessor.Gain = 0.001f;
- // ---- Reverb ----
+ intakeSoundProcessor = new SoundProcessor(sampleRate, intRadius * 2);
+ intakeSoundProcessor.Gain = 0.001f;
+
+ // Reverb
reverb = new OutdoorExhaustReverb(sampleRate);
- // Church: vast, highly reflective, bright
- reverb.DryMix = 1.0f; // always full dry signal
- reverb.EarlyMix = 0.5f; // distinct early reflections from distant walls
- reverb.TailMix = 0.9f; // huge tail, almost as loud as the dry sound
- reverb.Feedback = 0.9f; // long decay – roughly 3 s reverb time (with current delay lengths)
- reverb.DampingFreq = 6000f; // bright: high‑frequency energy stays for a long time
- reverb.MatrixCoeff = 0.5f; // default orthogonal mix
+ reverb.DryMix = 1.0f;
+ reverb.EarlyMix = 0.5f;
+ reverb.TailMix = 0.9f;
+ reverb.Feedback = 0.9f;
+ reverb.DampingFreq = 6000f;
- Console.WriteLine("=== EngineScenario (Stable) ===");
- Console.WriteLine($"Crankshaft inertia: {crankshaft.Inertia}");
- Console.WriteLine($"Pipe: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz");
+ Console.WriteLine("=== Engine with intake & cycle‑aware valves ===");
}
public override float Process()
{
- // ---- RPM governor: adjust throttle to maintain idle when no user input ----
- double currentRPM = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI);
- double throttle = Math.Clamp(Throttle, 0.05, 1.0); // never let it drop below a tiny value
-
- // ---- Target combustion pressure ----
+ double throttle = Math.Clamp(Throttle, 0.2, 1.0);
double targetPressure = throttle * FullLoadPeakPressure;
engineCyl.TargetPeakPressure = targetPressure;
- // ---- Simulate one timestep ----
engineCyl.Step(dt);
crankshaft.Step(dt);
- coupling.OrificeArea = engineCyl.OrificeArea;
+
+ couplingExhaust.OrificeArea = engineCyl.ExhaustOrificeArea;
+ couplingIntake.OrificeArea = engineCyl.IntakeOrificeArea;
+
solver.Step();
- // ---- Update exit port with safety clamps ----
- UpdateExitPort();
+ UpdateExhaustPort();
+ UpdateIntakePort();
+ float dryExhaust = exhaustSoundProcessor.Process(exhaustPort);
+ float dryIntake = intakeSoundProcessor.Process(intakePort);
+ float dry = dryExhaust + dryIntake;
- // ---- Generate audio ----
- float dry = soundProcessor.Process(exitPort);
float wet = reverb.Process(dry);
- time += dt;
- stepCount++;
+ if (++stepCount % LogInterval == 0) Log();
return wet;
}
- private void UpdateExitPort()
+ private void Log()
+ {
+ double rpm = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI);
+ double cycleDeg = (engineCyl.CycleAngle * 180.0 / Math.PI) % 720.0;
+ string stroke = cycleDeg < 180.0 ? "Power" :
+ cycleDeg < 360.0 ? "Exhaust" :
+ cycleDeg < 540.0 ? "Intake" : "Compression";
+
+ // Cylinder
+ double pCyl = engineCyl.Cylinder.Pressure;
+ double TCyl = engineCyl.Cylinder.Temperature;
+ double VCyl = engineCyl.Cylinder.Volume;
+ double mCyl = engineCyl.Cylinder.Mass;
+ double exhArea = engineCyl.ExhaustOrificeArea * 1e6; // mm²
+ double intArea = engineCyl.IntakeOrificeArea * 1e6; // mm²
+
+ // Exhaust pipe
+ int exhLast = exhaustPipe.GetCellCount() - 1;
+ double pExhEnd = exhaustPipe.GetCellPressure(exhLast);
+ double mdotExhOut = exhaustPipe.GetOpenEndMassFlow(); // positive out
+
+ // Intake pipe
+ double mdotIntIn = couplingIntake.LastMassFlowIntoVolume;
+ double pIntAmbEnd = intakePort.Pressure;
+
+ Console.WriteLine(
+ $"{stepCount,8} {stroke,-11} {cycleDeg,6:F1}° " +
+ $"RPM:{rpm,5:F0} " +
+ $"Cyl: p={pCyl/1e5,6:F3}bar T={TCyl,6:F0}K V={VCyl*1e6,6:F0}cm³ m={mCyl*1e3,6:F6}g " +
+ $"Valves: Exh={exhArea,5:F0}mm² Int={intArea,5:F0}mm² " +
+ $"Intake: p_end={pIntAmbEnd/1e5,6:F3}bar mdot_in={mdotIntIn,7:F4}kg/s " +
+ $"Exhaust: p_end={pExhEnd/1e5,6:F3}bar mdot_out={mdotExhOut,7:F4}kg/s");
+ }
+
+ private void UpdateExhaustPort()
{
int last = exhaustPipe.GetCellCount() - 1;
double p = exhaustPipe.GetCellPressure(last);
double rho = exhaustPipe.GetCellDensity(last);
double vel = exhaustPipe.GetCellVelocity(last);
- // Clamp density to physically possible values
- if (rho < 0.01) rho = 0.01; // never let it approach zero
- if (rho > 50.0) rho = 50.0; // never let it become absurd
-
- // Clamp velocity to ± 500 m/s (safe subsonic)
+ // Safety clamps
+ rho = Math.Clamp(rho, 0.01, 50.0);
vel = Math.Clamp(vel, -500.0, 500.0);
+ p = Math.Clamp(p, 1e4, 2e6);
- double outflowMassFlow = rho * vel * pipeArea;
+ double outflowMassFlow = rho * vel * exhPipeArea;
- // Clamp exit pressure to sensible range (0.1 – 20 bar)
- p = Math.Clamp(p, 1e4, 2e6);
-
- exitPort.Pressure = p;
- exitPort.Density = rho;
- exitPort.Temperature = p / (rho * 287.05);
- exitPort.MassFlowRate = -outflowMassFlow;
- exitPort.SpecificEnthalpy = 0.0;
+ exhaustPort.Pressure = p;
+ exhaustPort.Density = rho;
+ exhaustPort.Temperature = p / (rho * 287.05);
+ exhaustPort.MassFlowRate = -outflowMassFlow;
+ exhaustPort.SpecificEnthalpy = 0.0;
}
+ private void UpdateIntakePort()
+ {
+ // Use the actual valve mass flow (positive = into cylinder)
+ double mdotIntoEngine = couplingIntake.LastMassFlowIntoVolume;
+
+ // Use cylinder pressure/density for the port state (or intake pipe last cell)
+ double pCyl = engineCyl.Cylinder.Pressure;
+ double rhoCyl = engineCyl.Cylinder.Density;
+
+ intakePort.Pressure = Math.Max(pCyl, 100);
+ intakePort.Density = Math.Max(rhoCyl, 1e-6);
+ intakePort.Temperature = engineCyl.Cylinder.Temperature;
+ intakePort.MassFlowRate = mdotIntoEngine;
+ intakePort.SpecificEnthalpy = 0.0;
+ }
+
+ // ==================== Drawing ====================
public override void Draw(RenderWindow target)
{
float winW = target.GetView().Size.X;
@@ -169,10 +227,10 @@ namespace FluidSim.Core
return (float)Math.Clamp(t, -1.0, 1.0);
}
+ // ---- Cylinder ----
float cylW = 80f, cylH = 150f;
var cylRect = new RectangleShape(new Vector2f(cylW, cylH));
- cylRect.Position = new Vector2f(40f, centerY - cylH / 2f);
-
+ cylRect.Position = new Vector2f(200f, centerY - cylH / 2f);
double tempCyl = engineCyl.Cylinder.Temperature;
float tnCyl = NormaliseTemperature(tempCyl);
byte rC = (byte)(tnCyl > 0 ? 255 * tnCyl : 0);
@@ -181,33 +239,60 @@ namespace FluidSim.Core
cylRect.FillColor = new Color(rC, gC, bC);
target.Draw(cylRect);
- int n = exhaustPipe.GetCellCount();
- float pipeStartX = 120f, pipeEndX = winW - 60f;
- float pipeLen = pipeEndX - pipeStartX;
- float dx = pipeLen / (n - 1);
- float baseRadius = 20f;
+ // ---- Piston ----
+ float pistonWidth = cylW - 12f;
+ float pistonHeight = 16f;
+ float pistonFraction = (float)engineCyl.PistonPositionFraction;
+ float pistonTopY = cylRect.Position.Y + pistonFraction * (cylH - pistonHeight);
+ var pistonRect = new RectangleShape(new Vector2f(pistonWidth, pistonHeight))
+ {
+ Position = new Vector2f(cylRect.Position.X + 6f, pistonTopY),
+ FillColor = new Color(80, 80, 80)
+ };
+ target.Draw(pistonRect);
+
+ // ---- Exhaust pipe (rightwards) ----
+ DrawPipe(target, exhaustPipe, startX: 280f, endX: winW - 60f, centerY,
+ T_ambient, T_hot, T_cold, R, NormaliseTemperature, true);
+
+ // ---- Intake pipe (leftwards) ----
+ DrawPipe(target, intakePipe, startX: 200f, endX: 20f, centerY,
+ T_ambient, T_hot, T_cold, R, NormaliseTemperature, false);
+ }
+
+ private void DrawPipe(RenderWindow target, Pipe1D pipe,
+ float startX, float endX, float centerY,
+ float T_ambient, float T_hot, float T_cold, float R,
+ Func normaliseTemp, bool leftToRight)
+ {
+ int n = pipe.GetCellCount();
+ float dir = leftToRight ? 1f : -1f;
+ float pipeLen = Math.Abs(endX - startX);
+ float dx = pipeLen / (n - 1) * dir;
+ float baseRadius = leftToRight ? 20f : 14f; // exhaust thicker, intake thinner
var vertices = new Vertex[n * 2];
float ambPress = 101325f;
for (int i = 0; i < n; i++)
{
- float x = pipeStartX + i * dx;
- double p = exhaustPipe.GetCellPressure(i);
- double rho = exhaustPipe.GetCellDensity(i);
+ float x = startX + i * dx;
+ double p = pipe.GetCellPressure(i);
+ double rho = pipe.GetCellDensity(i);
double T = p / (rho * R);
float r = baseRadius * 0.3f * (float)(1.0 + (p - ambPress) / ambPress);
if (r < 2f) r = 2f;
- float tn = NormaliseTemperature(T);
- byte rCol = (byte)(tn > 0 ? 255 * tn : 0);
- byte bCol = (byte)(tn < 0 ? -255 * tn : 0);
- byte gCol = (byte)(255 * (1 - Math.Abs(tn)));
- var col = new Color(rCol, gCol, bCol);
+ float tn = normaliseTemp(T);
+ byte rC = (byte)(tn > 0 ? 255 * tn : 0);
+ byte bC = (byte)(tn < 0 ? -255 * tn : 0);
+ byte gC = (byte)(255 * (1 - Math.Abs(tn)));
+ var col = new Color(rC, gC, bC);
vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
}
+
target.Draw(vertices, PrimitiveType.TriangleStrip);
}
}