diff --git a/Components/Crankshaft.cs b/Components/Crankshaft.cs
new file mode 100644
index 0000000..99ea43e
--- /dev/null
+++ b/Components/Crankshaft.cs
@@ -0,0 +1,52 @@
+// Components/Crankshaft.cs
+using System;
+
+namespace FluidSim.Components
+{
+ public class Crankshaft
+ {
+ 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 Inertia { get; set; } = 0.2;
+ public double FrictionConstant { get; set; } = 2.0; // N·m
+ public double FrictionViscous { get; set; } = 0.005; // N·m per rad/s
+
+ 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;
+ CrankAngle = 0.0;
+ PreviousAngle = 0.0;
+ }
+
+ public void AddTorque(double torque) => externalTorque += torque;
+
+ public void Step(double dt)
+ {
+ // Save previous angle
+ 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
+
+ CrankAngle += AngularVelocity * dt;
+
+ // Wrap to [0, 4π)
+ if (CrankAngle >= 4.0 * Math.PI)
+ CrankAngle -= 4.0 * Math.PI;
+ else if (CrankAngle < 0)
+ CrankAngle += 4.0 * Math.PI;
+
+ externalTorque = 0.0;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/EngineCylinder.cs b/Components/EngineCylinder.cs
new file mode 100644
index 0000000..aa322fe
--- /dev/null
+++ b/Components/EngineCylinder.cs
@@ -0,0 +1,215 @@
+// EngineCylinder.cs (in Core namespace)
+using System;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+ public class EngineCylinder
+ {
+ public Volume0D Cylinder { get; private set; }
+ private Crankshaft crankshaft;
+
+ 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 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 burnDuration = 40.0 * Math.PI / 180.0;
+ private double targetBurnEnergy;
+ private double totalBurnMass;
+ private double preIgnitionMass, preIgnitionInternalEnergy;
+
+ private Random rand = new Random();
+ public double MisfireProbability { get; set; } = 0.02;
+ private bool misfireCurrent = false;
+
+ public int CombustionCount { get; private set; }
+ public int MisfireCount { get; private set; }
+
+ public EngineCylinder(Crankshaft crankshaft,
+ double bore, double stroke, double compressionRatio,
+ double pipeArea, int sampleRate)
+ {
+ this.crankshaft = crankshaft;
+ this.bore = bore;
+ this.stroke = stroke;
+ conRodLength = 2.0 * stroke;
+ this.compressionRatio = compressionRatio;
+ maxOrificeArea = pipeArea;
+ 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;
+ 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);
+
+ Cylinder = new Volume0D(V_clear, p_tdc, T_bdc * Math.Pow(V_bdc / V_clear, 1.4 - 1.0), sampleRate)
+ {
+ Gamma = 1.4,
+ GasConstant = 287.0
+ };
+ Cylinder.Volume = V_clear;
+ Cylinder.Mass = freshMass;
+ Cylinder.InternalEnergy = p_tdc * V_clear / (1.4 - 1.0);
+
+ preIgnitionMass = Cylinder.Mass;
+ preIgnitionInternalEnergy = Cylinder.InternalEnergy;
+ }
+
+ // ---- Piston kinematics (uses full cycle angle for position) ----
+ 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);
+ double sinT = Math.Sin(theta);
+ double L = conRodLength;
+
+ double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT);
+ double V = V_clear + pistonArea * s;
+
+ double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
+ double dVdθ = pistonArea * (R * sinT + (R * R * sinT * cosT) / sqrtTerm);
+ double dvdt = dVdθ * crankshaft.AngularVelocity;
+ return (V, dvdt);
+ }
+
+ // ---- Valve lift ----
+ private double ValveLift()
+ {
+ 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 rampFrac = ramp / duration;
+
+ if (t < rampFrac)
+ return t / rampFrac;
+ else if (t > 1.0 - rampFrac)
+ return (1.0 - t) / rampFrac;
+ else
+ return 1.0;
+ }
+
+ private double WiebeFraction(double angleFromIgnition)
+ {
+ if (angleFromIgnition >= burnDuration) return 1.0;
+ double a = 5.0, m = 2.0;
+ double x = angleFromIgnition / burnDuration;
+ return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1));
+ }
+
+ // ---- Torque from pressure ----
+ private double ComputeTorque()
+ {
+ double p = Cylinder.Pressure;
+ double ambient = 101325.0;
+ double force = (p - ambient) * pistonArea;
+ if (force <= 0) return 0.0;
+
+ double theta = crankshaft.CrankAngle % (2.0 * Math.PI);
+ double R = stroke / 2.0;
+ double L = conRodLength;
+ double sinT = Math.Sin(theta);
+ double cosT = Math.Cos(theta);
+
+ double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
+ double lever = R * sinT * (1.0 + (R * cosT) / sqrtTerm);
+ return force * lever;
+ }
+
+ 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)
+
+ 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);
+
+ preIgnitionMass = Cylinder.Mass;
+ preIgnitionInternalEnergy = Cylinder.InternalEnergy;
+
+ if (misfireCurrent)
+ {
+ MisfireCount++;
+ }
+ else
+ {
+ double V = V_clear;
+ targetBurnEnergy = TargetPeakPressure * V / (gamma - 1.0);
+ totalBurnMass = TargetPeakPressure * V / (287.0 * PeakTemperature);
+ burnInProgress = true;
+ burnStartAngle = cycleAngle;
+ CombustionCount++;
+ }
+ }
+
+ // ----- Burn progress -----
+ if (burnInProgress)
+ {
+ double angleFromIgnition = cycleAngle - burnStartAngle;
+ if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI; // wrap if needed
+
+ 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;
+ }
+ }
+
+ // ----- Piston motion -----
+ var (vol, dvdt) = PistonKinematics(cycleAngle);
+ Cylinder.Volume = vol;
+ Cylinder.Dvdt = dvdt;
+
+ // ----- Torque contribution -----
+ double torque = ComputeTorque();
+ crankshaft.AddTorque(torque);
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
index ac444a7..6c414cb 100644
--- a/Components/Pipe1D.cs
+++ b/Components/Pipe1D.cs
@@ -53,6 +53,10 @@ namespace FluidSim.Components
// Pre‑computed for damping
private float _laminarCoeff; // 8*mu / r^2, then multiplied by DampingMultiplier
+ // ---- Energy loss (Newton cooling) ----
+ private float _ambientEnergyReference; // total energy density at ambient (Pamb / (γ-1))
+ public float EnergyRelaxationRate { get; set; } = 0.0f; // 1/s
+
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
{
float dtGlobal = 1f / sampleRate;
@@ -89,12 +93,14 @@ namespace FluidSim.Components
float radius = _diameter * 0.5f;
_laminarCoeff = 8f * mu_air / (radius * radius); // will be multiplied by DampingMultiplier at each step
+ // Ambient reference energy (internal energy per unit volume at 101325 Pa)
+ _ambientEnergyReference = 101325f / (_gamma - 1f); // ≈ 253312.5 J/m³
+
PortA = new Port();
PortB = new Port();
}
- // ==================== PUBLIC API (unchanged) ============================
-
+ // ==================== PUBLIC API ============================
public void SetABoundaryType(BoundaryType type) => _aBCType = type;
public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
@@ -237,8 +243,6 @@ namespace FluidSim.Components
// --- 2. Internal faces (1 .. n-1) – SIMD ---------------------------
int vectorSize = Vector.Count;
int lastSimdFace = n - vectorSize; // highest face index that starts a full vector block
- // Face index f is between cell f-1 (left) and cell f (right).
- // We want to cover faces 1..n-1.
for (int f = 1; f <= lastSimdFace; f += vectorSize)
{
SimdInternalFluxBlock(f, vectorSize);
@@ -262,7 +266,7 @@ namespace FluidSim.Components
float pR = PressureScalar(n - 1);
ComputeRightBoundaryFlux(rhoR, uR, pR, out _fluxM[n], out _fluxP[n], out _fluxE[n]);
- // --- 4. Cell update + damping – SIMD ------------------------------
+ // --- 4. Cell update + damping + energy loss – SIMD -----------------
SimdCellUpdate(dt);
}
@@ -417,14 +421,9 @@ namespace FluidSim.Components
// ==================== SIMD INTERNAL FACE ROUTINE ========================
private void SimdInternalFluxBlock(int startFace, int count)
{
- // startFace is the first face index; we process 'count' consecutive faces.
- // For face f, left cell = f-1, right cell = f.
- // We load left and right states for faces [startFace .. startFace+count-1].
-
int leftIdx = startFace - 1;
int rightIdx = startFace;
- // Load conserved variables for left cells and right cells as vectors.
Vector rL = new Vector(_rho, leftIdx);
Vector ruL = new Vector(_rhou, leftIdx);
Vector EL = new Vector(_E, leftIdx);
@@ -433,7 +432,6 @@ namespace FluidSim.Components
Vector ruR = new Vector(_rhou, rightIdx);
Vector ER = new Vector(_E, rightIdx);
- // Derived quantities: u = ru / r, p = (gamma-1)*(E - 0.5*ru^2 / r)
Vector uL = ruL / rL;
Vector uR = ruR / rR;
@@ -444,35 +442,30 @@ namespace FluidSim.Components
Vector pL = gammaMinus1 * (EL - half * ruL * ruL / rL);
Vector pR = gammaMinus1 * (ER - half * ruR * ruR / rR);
- // Sound speeds
Vector cL = Vector.SquareRoot(gammaVec * pL / rL);
Vector cR = Vector.SquareRoot(gammaVec * pR / rR);
- // Wave speeds
Vector SL = Vector.Min(uL - cL, uR - cR);
Vector SR = Vector.Max(uL + cL, uR + cR);
- // Star speed
Vector num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR);
Vector den = rL * (SL - uL) - rR * (SR - uR);
Vector Ss = num / den;
- // Total energy per unit mass (E/rho) for left/right (needed for star region)
Vector eL = EL / rL;
Vector eR = ER / rR;
- // --- Compute all four possible flux vectors ---
// Left flux
Vector Fm_L = ruL;
Vector Fp_L = ruL * uL + pL;
- Vector Fe_L = (EL + pL) * uL; // (r*E + p)*u
+ Vector Fe_L = (EL + pL) * uL;
// Right flux
Vector Fm_R = ruR;
Vector Fp_R = ruR * uR + pR;
Vector Fe_R = (ER + pR) * uR;
- // Star‑left fluxes (when SL < 0 < Ss)
+ // Star‑left fluxes
Vector diffL = SL - uL;
Vector dL_den = SL - Ss;
Vector rsL = rL * diffL / dL_den;
@@ -482,7 +475,7 @@ namespace FluidSim.Components
Vector Fp_starL = rsL * Ss * Ss + psSL;
Vector Fe_starL = (rsL * EsL + psSL) * Ss;
- // Star‑right fluxes (when Ss < 0 < SR)
+ // Star‑right fluxes
Vector diffR = SR - uR;
Vector dR_den = SR - Ss;
Vector rsR = rR * diffR / dR_den;
@@ -492,14 +485,12 @@ namespace FluidSim.Components
Vector Fp_starR = rsR * Ss * Ss + psSR;
Vector Fe_starR = (rsR * EsR + psSR) * Ss;
- // --- Select the correct flux based on wave signs ---
var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector.Zero);
var maskSRle0 = Vector.LessThanOrEqual(SR, Vector.Zero);
- var maskMiddle = ~(maskSLge0 | maskSRle0); // SL<0 && SR>0
+ var maskMiddle = ~(maskSLge0 | maskSRle0);
var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector.Zero);
var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector.Zero);
- // Start with left flux, override with right/star as needed
Vector fm = Vector.ConditionalSelect(maskSLge0, Fm_L,
Vector.ConditionalSelect(maskSRle0, Fm_R,
Vector.ConditionalSelect(maskStarL, Fm_starL,
@@ -515,13 +506,12 @@ namespace FluidSim.Components
Vector.ConditionalSelect(maskStarL, Fe_starL,
Vector.ConditionalSelect(maskStarR, Fe_starR, Vector.Zero))));
- // Store to flux arrays at indices startFace .. startFace+count-1
fm.CopyTo(_fluxM, startFace);
fp.CopyTo(_fluxP, startFace);
fe.CopyTo(_fluxE, startFace);
}
- // ==================== SIMD CELL UPDATE + DAMPING ========================
+ // ==================== SIMD CELL UPDATE + DAMPING + ENERGY LOSS =========
private void SimdCellUpdate(float dt)
{
float dt_dx = dt / _dx;
@@ -534,9 +524,11 @@ namespace FluidSim.Components
int n = _n;
int lastSimdCell = n - vectorSize;
- // Pre‑define constants used in clamping
+ // Pre‑defined constants used in clamping
Vector half = new Vector(0.5f);
Vector gammaMinus1 = new Vector(_gamma - 1f);
+ Vector ambientEnergyVec = new Vector(_ambientEnergyReference);
+ Vector energyRelaxRateVec = new Vector(EnergyRelaxationRate);
for (int i = 0; i <= lastSimdCell; i += vectorSize)
{
@@ -559,8 +551,12 @@ namespace FluidSim.Components
Vector newE = E - vDtDx * (flxE_R - flxE_L);
// Damping
- Vector factor = Vector.Exp(-vCoeff / r * vDt);
- newRu *= factor;
+ Vector dampingFactor = Vector.Exp(-vCoeff / r * vDt);
+ newRu *= dampingFactor;
+
+ // Energy loss (Newton cooling toward ambient)
+ Vector relaxFactor = Vector.Exp(-energyRelaxRateVec * vDt);
+ newE = ambientEnergyVec + (newE - ambientEnergyVec) * relaxFactor;
// Clamp density
newR = Vector.Max(newR, new Vector(1e-12f));
@@ -576,6 +572,7 @@ namespace FluidSim.Components
}
// Scalar remainder
+ float relaxRate = EnergyRelaxationRate;
for (int i = Math.Max(0, lastSimdCell + 1); i < n; i++)
{
float r = _rho[i];
@@ -584,15 +581,19 @@ namespace FluidSim.Components
float dM = _fluxM[i + 1] - _fluxM[i];
float dP = _fluxP[i + 1] - _fluxP[i];
- float dE = _fluxE[i + 1] - _fluxE[i];
+ float dE_flux = _fluxE[i + 1] - _fluxE[i];
float newR = r - dt_dx * dM;
float newRu = ru - dt_dx * dP;
- float newE = E - dt_dx * dE;
+ float newE = E - dt_dx * dE_flux;
// Damping
- float factor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
- newRu *= factor;
+ float dampingFactor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
+ newRu *= dampingFactor;
+
+ // Energy loss
+ float relaxFactor = MathF.Exp(-relaxRate * dt);
+ newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
// Clamps
if (newR < 1e-12f) newR = 1e-12f;
diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs
index 9c48efb..3fb99bd 100644
--- a/Core/SoundProcessor.cs
+++ b/Core/SoundProcessor.cs
@@ -4,14 +4,20 @@ namespace FluidSim.Core
{
public class SoundProcessor
{
- // Monopole state
- private double lastMassFlow = 0.0;
private double dt;
+ private double pipeArea;
+ private double ambientPressure = 101325.0;
- // Resonant band‑pass filter (second‑order)
- private double b0, b1, b2, a1, a2;
- private double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
- private double pipeLength;
+ // Monopole source state
+ private double lastMassFlow = 0.0;
+
+ // Gains
+ private float masterGain = 0.0005f;
+ private float pressureGain = 0.12f;
+ private float turbulenceGain = 0.05f;
+ private float turbulence = 0.05f;
+
+ private PinkNoiseGenerator pinkNoise;
// Reverb (outdoor)
private float[] delayLine;
@@ -20,36 +26,11 @@ namespace FluidSim.Core
private float lowpassCoeff = 0.70f;
private float lastFeedbackSample = 0f;
- // Turbulence (pink noise scaled by U³)
- private PinkNoiseGenerator pinkNoise;
- private float turbulenceGain = 0.05f;
- private double pipeArea;
- private double ambientPressure = 101325.0;
-
- // Gains
- private float masterGain = 0.0005f;
- private float pressureGain = 0.12f;
-
- public SoundProcessor(int sampleRate, double pipeLengthMeters, double pipeDiameterMeters = 0.04, float reverbTimeMs = 200.0f)
+ public SoundProcessor(int sampleRate, double pipeDiameterMeters, float reverbTimeMs = 200.0f)
{
dt = 1.0 / sampleRate;
- pipeLength = pipeLengthMeters;
pipeArea = Math.PI * Math.Pow(pipeDiameterMeters / 2.0, 2.0);
- // Design resonant filter at pipe fundamental frequency
- double c = 340.0;
- double f0 = c / (4.0 * pipeLength);
- double Q = 15.0;
- double omega = 2.0 * Math.PI * f0;
- double alpha = Math.Sin(omega * dt) / (2.0 * Q);
- double norm = 1.0 / (1.0 + alpha);
- b0 = alpha * norm;
- b1 = 0.0;
- b2 = -alpha * norm;
- a1 = -2.0 * Math.Cos(omega * dt) * norm;
- a2 = (1.0 - alpha) * norm;
-
- // Reverb delay line
int delaySamples = (int)(sampleRate * reverbTimeMs / 1000.0);
delayLine = new float[delaySamples];
writeIndex = 0;
@@ -62,65 +43,60 @@ namespace FluidSim.Core
get => masterGain;
set => masterGain = value;
}
-
public float PressureGain
{
get => pressureGain;
set => pressureGain = value;
}
-
public float TurbulenceGain
{
get => turbulenceGain;
set => turbulenceGain = value;
}
+ public float Turbulence
+ {
+ get => turbulence;
+ set => turbulence = value;
+ }
public void SetAmbientPressure(double p) => ambientPressure = p;
- public void SetPipeDiameter(double diameterMeters) => pipeArea = Math.PI * Math.Pow(diameterMeters / 2.0, 2.0);
+ public void SetPipeDiameter(double diameterMeters) =>
+ pipeArea = Math.PI * Math.Pow(diameterMeters / 2.0, 2.0);
public float Process(float massFlow, float pipeEndPressure)
{
- // 1. Monopole source: d(mdot)/dt
+ // 1. Monopole: d(mdot)/dt
double derivative = (massFlow - lastMassFlow) / dt;
- derivative = Math.Clamp(derivative, -500, 500);
lastMassFlow = massFlow;
float monopole = (float)(derivative * masterGain);
- // 2. Pressure difference (low‑frequency component)
+ // 2. Pressure component
float pressureDiff = (float)((pipeEndPressure - ambientPressure) / ambientPressure) * pressureGain;
+
float mixed = monopole + pressureDiff;
- // DO NOT clamp here – let the filter and final clamp handle dynamics
- // 3. Resonant band‑pass filter
- double y = b0 * mixed + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;
- x2 = x1; x1 = mixed;
- y2 = y1; y1 = y;
- float resonant = (float)Math.Clamp(y, -1f, 1f);
-
- // 4. Turbulence noise: amplitude ∝ U³ (empirical for low speeds)
+ // 3. Turbulence: amplitude ∝ U^8
double velocity = massFlow / (pipeArea * 1.225);
- double Uref = 100.0;
- double turbulenceAmp = Math.Pow(Math.Abs(velocity) / Uref, 3.0);
+ double turbulenceAmp = Math.Pow(Math.Abs(velocity) * turbulence, 3.0);
float pink = pinkNoise.Next() * turbulenceGain * (float)turbulenceAmp;
- resonant += pink;
- resonant = Math.Clamp(resonant, -1f, 1f);
+ float combined = mixed + pink;
- // 5. Outdoor reverb
+ // 4. Outdoor reverb
float delayed = delayLine[writeIndex];
float filteredDelay = delayed * lowpassCoeff + lastFeedbackSample * (1f - lowpassCoeff);
lastFeedbackSample = filteredDelay;
float wet = delayed + filteredDelay * feedback;
- delayLine[writeIndex] = resonant + filteredDelay * feedback;
+ delayLine[writeIndex] = combined + filteredDelay * feedback;
writeIndex = (writeIndex + 1) % delayLine.Length;
- // 6. Dry/wet mix
- float output = resonant * 0.7f + wet * 0.3f;
- output = MathF.Tanh(output);
- return output;
+ // 5. Dry/wet mix
+ float output = combined * 0.7f + wet * 0.3f;
+ return MathF.Tanh(output);
}
}
+ // PinkNoiseGenerator unchanged, same as before
internal class PinkNoiseGenerator
{
private readonly Random random = new Random();
diff --git a/Program.cs b/Program.cs
index 11a2048..201c0cc 100644
--- a/Program.cs
+++ b/Program.cs
@@ -12,24 +12,27 @@ public class Program
private const double DrawFrequency = 60.0;
private static Scenario scenario;
- // Speed control
- private static double desiredSpeed = 0.0001;
- //private static double desiredSpeed = 1;
+ // Speed control (existing + new throttle)
+ private static double desiredSpeed = 0.01;
private static double currentSpeed = desiredSpeed;
private const double MinSpeed = 0.0001;
private const double MaxSpeed = 1.0;
private const double ScrollFactor = 1.1;
- // Space‑toggle state
- private static double lastDesiredSpeed = 0.1; // remembers the last non‑1.0 speed
- private static bool isRealTime = false; // starts in slow‑mo (desiredSpeed != 1.0)
+ private static double lastDesiredSpeed = 0.1;
+ private static bool isRealTime = false;
+
+ // 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 static volatile bool running = true;
public static void Main()
{
var mode = new VideoMode(new Vector2u(1280, 720));
- var window = new RenderWindow(mode, "FluidSim - Helmholtz Resonator");
+ var window = new RenderWindow(mode, "FluidSim - Engine (W = throttle)");
window.SetVerticalSyncEnabled(true);
window.Closed += (_, _) => { running = false; window.Close(); };
window.MouseWheelScrolled += OnMouseWheel;
@@ -39,12 +42,7 @@ public class Program
soundEngine.Volume = 70;
soundEngine.Start();
- // Choose one scenario. The Helmholtz resonator is fully updated.
- //scenario = new HelmholtzResonatorScenario();
- //scenario = new PipeResonatorScenario(); // needs update to new API
- //scenario = new SodShockTubeScenario(); // needs update to new API
- scenario = new EngineScenario(); // also works (provided earlier)
-
+ scenario = new EngineScenario();
scenario.Initialize(SampleRate);
var stopwatch = Stopwatch.StartNew();
@@ -52,18 +50,13 @@ public class Program
double drawInterval = 1.0 / DrawFrequency;
double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
- // Resampling buffer
var simBuffer = new List(4096);
double readIndex = 0.0;
-
- // Prime the buffer with a few samples
for (int i = 0; i < 4; i++)
simBuffer.Add(scenario.Process());
long totalSimSteps = simBuffer.Count;
long totalOutputSamples = 0;
-
- double lastRealTime = stopwatch.Elapsed.TotalSeconds;
const int outputChunk = 256;
float[] outputBuf = new float[outputChunk];
@@ -72,16 +65,22 @@ public class Program
window.DispatchEvents();
double currentRealTime = stopwatch.Elapsed.TotalSeconds;
- double dtSpeed = currentRealTime - lastSpeedUpdateTime;
+ double dtClock = currentRealTime - lastSpeedUpdateTime;
lastSpeedUpdateTime = currentRealTime;
- // Smoothly transition currentSpeed → desiredSpeed
- double smoothingRate = 8.0; // higher = faster catch‑up
- currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
+ // Smooth simulation speed
+ double speedSmoothing = 8.0;
+ currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-speedSmoothing * dtClock));
- // ---------- Generate audio ----------
+ // ---- THROTTLE INPUT ----
+ targetThrottle = Keyboard.IsKeyPressed(Keyboard.Key.W) ? 1.0 : 0.0;
+ currentThrottle += (targetThrottle - currentThrottle) * (1.0 - Math.Exp(-ThrottleSmoothing * dtClock));
+ // Push to engine scenario (if it's an EngineScenario)
+ if (scenario is EngineScenario engine)
+ engine.Throttle = currentThrottle;
+
+ // Generate audio
double targetAudioClock = currentRealTime + 0.05;
-
while (totalOutputSamples < targetAudioClock * SampleRate && running)
{
int toGenerate = (int)Math.Min(
@@ -103,11 +102,9 @@ public class Program
int i0 = (int)readIndex;
int i1 = i0 + 1;
double frac = readIndex - i0;
-
float y0 = simBuffer[Math.Clamp(i0, 0, simBuffer.Count - 1)];
float y1 = simBuffer[Math.Clamp(i1, 0, simBuffer.Count - 1)];
outputBuf[i] = (float)(y0 + (y1 - y0) * frac);
-
readIndex += currentSpeed;
while (readIndex >= 1.0 && simBuffer.Count > 2)
@@ -119,23 +116,23 @@ public class Program
int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
totalOutputSamples += accepted;
-
if (accepted < toGenerate)
break;
}
- // ---------- Drawing & window title ----------
+ // Drawing & title
if (currentRealTime - lastDrawTime >= drawInterval)
{
double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate);
double simTime = totalSimSteps / (double)SampleRate;
string toggleHint = isRealTime ? "[Space] slow mo" : "[Space] real time";
+ string throttleHint = Keyboard.IsKeyPressed(Keyboard.Key.W) ? "W held" : "W released";
window.SetTitle(
- $"{toggleHint} Sim: {simTime:F2}s | " +
- $"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
- $"Actual: {actualSpeed:F2}x"
+ $"{toggleHint} {throttleHint} " +
+ $"Thr: {currentThrottle:F2} " +
+ $"Speed: {currentSpeed:F3}x → {desiredSpeed:F3}x " +
+ $"Act: {actualSpeed:F2}x"
);
-
window.Clear(Color.Black);
scenario.Draw(window);
window.Display();
@@ -147,22 +144,18 @@ public class Program
window.Dispose();
}
+ // (Mouse wheel, space toggle unchanged)
private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
{
bool wasRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-
if (e.Delta > 0)
desiredSpeed *= ScrollFactor;
else if (e.Delta < 0)
desiredSpeed /= ScrollFactor;
desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed);
-
- // Update the remembered slow‑mo speed (unless we are exactly at 1.0)
if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6)
lastDesiredSpeed = desiredSpeed;
-
- // Update isRealTime flag
isRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
}
@@ -171,15 +164,9 @@ public class Program
if (e.Code == Keyboard.Key.Space)
{
if (isRealTime)
- {
- // Switch to the remembered slow speed
desiredSpeed = lastDesiredSpeed;
- }
else
- {
- // Switch back to real time
desiredSpeed = 1.0;
- }
isRealTime = !isRealTime;
}
}
diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs
index 01cd681..769ac2e 100644
--- a/Scenarios/EngineScenario.cs
+++ b/Scenarios/EngineScenario.cs
@@ -8,7 +8,8 @@ namespace FluidSim.Core
public class EngineScenario : Scenario
{
private Solver solver;
- private Volume0D cylinder;
+ private Crankshaft crankshaft;
+ private EngineCylinder engineCyl;
private Pipe1D exhaustPipe;
private PipeVolumeConnection coupling;
private SoundProcessor soundProcessor;
@@ -16,259 +17,98 @@ namespace FluidSim.Core
private double dt;
private double ambientPressure = 101325.0;
private double time;
-
- // ---- 4‑stroke cycle angle (0 … 4π) ----
- private double cycleCrankAngle = 0.0; // 0 to 4π, then resets
- private const double TargetRPM = 1000.0;
- private double angularVelocity; // rad/s of crankshaft
-
- // ---- Engine geometry ----
- private double bore = 0.065; // 65 mm
- private double stroke = 0.0565; // 56.5 mm → 250 cc
- private double conRodLength = 0.113; // roughly 2 * stroke
- private double compressionRatio = 10.0;
- private double V_disp; // displacement volume
- private double V_clear; // clearance volume
-
- // ---- Combustion ----
- private const double CombustionPressure = 50.0 * 101325.0;
- private const double CombustionTemperature = 2500.0;
- private bool burnInProgress = false;
- private double burnStartAngle; // cycle angle when ignition began
- private const double BurnDurationDeg = 40.0;
- private const double BurnDurationRad = BurnDurationDeg * Math.PI / 180.0;
- private double targetBurnEnergy;
- private double totalBurnMass;
-
- // Pre‑ignition state (compressed fresh charge) for misfire restoration
- private double preIgnitionMass;
- private double preIgnitionInternalEnergy;
-
- // ---- Valve timing ----
- private const double ValveOpenStart = 120.0 * Math.PI / 180.0; // 120° after TDC power
- private const double ValveOpenEnd = 480.0 * Math.PI / 180.0; // 480° ≈ 120° after TDC exhaust
- private const double ValveRampWidth = 30.0 * Math.PI / 180.0; // 30° ramps
-
- private double maxOrificeArea;
-
- // ---- Misfire ----
- private Random rand = new Random();
- private const double MisfireProbability = 0.02;
- private bool isMisfiring = false;
-
- // ---- Logging ----
private int stepCount = 0;
- private const int LogStepInterval = 10000;
- private int combustionCount = 0;
- private int misfireCount = 0;
+ private const int LogInterval = 10000;
+
+ // Throttle 0..1 → target combustion pressure
+ public double Throttle { get; set; } = 0.05; // tiny throttle to keep idle
+ private const double IdlePeakPressure = 5.0 * 101325.0; // 5 bar
+ private const double MaxPeakPressure = 50.0 * 101325.0; // 50 bar
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
- angularVelocity = TargetRPM * 2.0 * Math.PI / 60.0;
- // Displacement volume
- V_disp = (Math.PI / 4.0) * bore * bore * stroke;
- V_clear = V_disp / (compressionRatio - 1.0);
-
- // Cylinder (starts at TDC clearance volume with compressed ambient charge)
- double initialPressure = ambientPressure * Math.Pow(compressionRatio, 1.4); // isentropic compression
- double initialTemperature = 300.0 * Math.Pow(compressionRatio, 1.4 - 1.0);
- double initialVolume = V_clear;
- cylinder = new Volume0D(initialVolume, initialPressure, initialTemperature, sampleRate)
+ // Crankshaft (inertia + friction)
+ crankshaft = new Crankshaft(initialRPM: 100.0) // starter speed
{
- Gamma = 1.4,
- GasConstant = 287.0
+ Inertia = 0.05,
+ FrictionConstant = 1.0,
+ FrictionViscous = 0.01
};
- // Exhaust pipe (2.5 m long, 3 cm radius)
- double pipeLength = 2.5;
- double pipeRadius = 0.03;
+ // Pipe
+ double pipeLength = 0.5;
+ double pipeRadius = 0.1;
double pipeArea = Math.PI * pipeRadius * pipeRadius;
- maxOrificeArea = pipeArea;
- exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 100);
+ exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 60);
exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure);
+ exhaustPipe.EnergyRelaxationRate = 0f;
+ exhaustPipe.DampingMultiplier = 0;
- // Coupling (valve initially closed)
- coupling = new PipeVolumeConnection(cylinder, exhaustPipe, true, orificeArea: 0.0);
+ // Cylinder (coupled to crankshaft)
+ engineCyl = new EngineCylinder(crankshaft,
+ bore: 0.065, stroke: 0.0565, compressionRatio: 10.0,
+ pipeArea: pipeArea, sampleRate: sampleRate);
+ // Coupling (valve → pipe)
+ coupling = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0);
+
+ // Solver
solver = new Solver();
solver.SetTimeStep(dt);
- solver.AddVolume(cylinder);
+ solver.AddVolume(engineCyl.Cylinder);
solver.AddPipe(exhaustPipe);
solver.AddConnection(coupling);
- // Open end with characteristic radiation (softer reflections)
solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, ambientPressure);
- // Sound processor (keep your carefully tuned gains)
- soundProcessor = new SoundProcessor(sampleRate, pipeLength, pipeRadius * 2, reverbTimeMs: 200.0f);
- soundProcessor.MasterGain = 0.0002f;
- soundProcessor.PressureGain = 10.0f;
- soundProcessor.TurbulenceGain = 0.00005f;
+ // Sound (your tuned gains)
+ soundProcessor = new SoundProcessor(sampleRate, pipeRadius * 2, reverbTimeMs: 500.0f);
+ soundProcessor.MasterGain = 0.0f; //0.00001f;
+ soundProcessor.PressureGain = 0.1f;
+ soundProcessor.TurbulenceGain = 0.0f;
+ soundProcessor.Turbulence = 0.001f;
soundProcessor.SetAmbientPressure(ambientPressure);
- // Log startup info
- Console.WriteLine("=== EngineScenario (improved physics) ===");
- Console.WriteLine($"Target RPM: {TargetRPM}, Misfire prob: {MisfireProbability * 100:F1}%");
- Console.WriteLine($"Bore x Stroke: {bore*1000:F0} x {stroke*1000:F0} mm, CR: {compressionRatio:F1}");
- Console.WriteLine($"Pipe length: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz");
- Console.WriteLine($"Combustion: {CombustionPressure/101325:F0} bar, {CombustionTemperature} K");
- Console.WriteLine($"Valve opens at {ValveOpenStart*180/Math.PI:F0}°, closes at {ValveOpenEnd*180/Math.PI:F0}° (ramp {ValveRampWidth*180/Math.PI:F0}°)");
- Console.WriteLine($"Burn duration: {BurnDurationDeg}°");
- Console.WriteLine("Time[s] Crank[°] Vol[cc] MassFlow[kg/s] Comb# Misfire");
- Console.WriteLine("-------------------------------------------------------------");
- }
-
- // ---- Piston volume & derivative ----
- private (double volume, double dvdt) PistonKinematics(double theta)
- {
- // theta = crankshaft angle (0 at TDC of power stroke)
- double R = stroke / 2.0;
- double cosT = Math.Cos(theta);
- double sinT = Math.Sin(theta);
- double L = conRodLength;
-
- // Slider‑crank position relative to TDC
- double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT);
- double V = V_clear + (Math.PI / 4.0) * bore * bore * s;
-
- // Derivative dV/dθ
- double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
- double dVdθ = (Math.PI / 4.0) * bore * bore * (R * sinT + (R * R * sinT * cosT) / sqrtTerm);
-
- double dvdt = dVdθ * angularVelocity; // rad/s → volume change rate
- return (V, dvdt);
- }
-
- // ---- Valve lift (trapezoidal) ----
- private double ValveOpenRatio(double cycleRad)
- {
- // cycleRad: 0 … 4π
- if (cycleRad < ValveOpenStart || cycleRad > ValveOpenEnd)
- return 0.0;
-
- double duration = ValveOpenEnd - ValveOpenStart;
- double ramp = ValveRampWidth;
-
- double t = (cycleRad - ValveOpenStart) / duration;
-
- if (t < ramp / duration)
- return t / (ramp / duration);
- else if (t > 1.0 - ramp / duration)
- return (1.0 - t) / (ramp / duration);
- else
- return 1.0;
- }
-
- // ---- Wiebe burn fraction ----
- private double WiebeFraction(double angleFromIgnition)
- {
- if (angleFromIgnition >= BurnDurationRad) return 1.0;
- double a = 5.0, m = 2.0;
- double x = angleFromIgnition / BurnDurationRad;
- return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1));
+ Console.WriteLine("=== EngineScenario (torque‑driven RPM, throttle = pressure) ===");
+ Console.WriteLine($"Crankshaft inertia: {crankshaft.Inertia}, friction: {crankshaft.FrictionConstant} + {crankshaft.FrictionViscous}*ω");
+ Console.WriteLine($"Throttle range: {IdlePeakPressure/101325:F0} – {MaxPeakPressure/101325:F0} bar");
+ Console.WriteLine($"Pipe: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz");
}
public override float Process()
{
- // Advance cycle crank angle
- cycleCrankAngle += angularVelocity * dt;
- if (cycleCrankAngle >= 4.0 * Math.PI) // 720°
- {
- cycleCrankAngle -= 4.0 * Math.PI;
- isMisfiring = rand.NextDouble() < MisfireProbability;
+ // 1. Map throttle to target peak pressure
+ double targetPressure = IdlePeakPressure + Throttle * (MaxPeakPressure - IdlePeakPressure);
+ engineCyl.TargetPeakPressure = targetPressure;
- // ---- Prepare cylinder for new power stroke ----
- // Fill cylinder with fresh charge at BDC, then compress isentropically to TDC clearance volume.
- double T_bdc = 300.0; // intake temperature
- double p_bdc = ambientPressure; // intake pressure
- double V_bdc = V_clear + V_disp; // volume at BDC (intake valve closing)
- double freshMass = p_bdc * V_bdc / (287.0 * T_bdc);
- double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0);
+ // 2. Step the cylinder (adds torque to crankshaft, updates valve)
+ engineCyl.Step(dt);
- // Compress isentropically to V_clear
- double V1 = V_bdc, V2 = V_clear;
- double gamma = 1.4;
- double p2 = p_bdc * Math.Pow(V1 / V2, gamma);
- double T2 = T_bdc * Math.Pow(V1 / V2, gamma - 1);
+ // 3. Integrate crankshaft (applies friction, updates RPM)
+ crankshaft.Step(dt);
- // Set compressed state
- cylinder.Volume = V_clear;
- cylinder.Mass = freshMass;
- cylinder.InternalEnergy = p2 * V_clear / (gamma - 1.0); // consistent with pressure/temperature
+ // 4. Set orifice area for coupling
+ coupling.OrificeArea = engineCyl.OrificeArea;
- // Store pre‑ignition state for misfire recovery
- preIgnitionMass = cylinder.Mass;
- preIgnitionInternalEnergy = cylinder.InternalEnergy;
-
- if (isMisfiring)
- {
- // No combustion – just expand from compressed state
- misfireCount++;
- }
- else
- {
- // Start Wiebe burn
- double V = V_clear;
- targetBurnEnergy = CombustionPressure * V / (gamma - 1.0);
- double R = 287.0;
- totalBurnMass = CombustionPressure * V / (R * CombustionTemperature);
- burnInProgress = true;
- burnStartAngle = cycleCrankAngle; // now = 0
- combustionCount++;
- }
- }
-
- // ---- Combustion progress (if active) ----
- if (burnInProgress)
- {
- double angleFromIgnition = cycleCrankAngle - burnStartAngle;
- if (angleFromIgnition >= BurnDurationRad)
- {
- // Burn complete
- cylinder.Mass = totalBurnMass;
- cylinder.InternalEnergy = targetBurnEnergy;
- burnInProgress = false;
- }
- else
- {
- double fraction = WiebeFraction(angleFromIgnition);
- // Interpolate between pre‑ignition (compressed charge) and final burned state
- double gamma = 1.4;
- double V = cylinder.Volume; // still near clearance
- double baseEnergy = preIgnitionInternalEnergy;
- double baseMass = preIgnitionMass;
-
- cylinder.InternalEnergy = baseEnergy * (1.0 - fraction) + targetBurnEnergy * fraction;
- cylinder.Mass = baseMass * (1.0 - fraction) + totalBurnMass * fraction;
- }
- }
-
- // ---- Update cylinder volume from piston kinematics ----
- double theta = cycleCrankAngle % (2.0 * Math.PI); // crank angle for piston position
- var (vol, dvdt) = PistonKinematics(theta);
- cylinder.Volume = vol;
- cylinder.Dvdt = dvdt;
-
- // ---- Valve lift & orifice area ----
- double lift = ValveOpenRatio(cycleCrankAngle);
- coupling.OrificeArea = maxOrificeArea * lift;
-
- // ---- Solver step ----
+ // 5. Fluid solver step
float massFlow = solver.Step();
float endPressure = (float)exhaustPipe.GetCellPressure(exhaustPipe.GetCellCount() - 1);
- // ---- Audio (no filter, feed raw pressure) ----
+ // 6. Audio
float audioSample = soundProcessor.Process(massFlow, endPressure);
- // Log occasionally
time += dt;
stepCount++;
- if (stepCount % LogStepInterval == 0)
+ if (stepCount % LogInterval == 0) {
+ Console.WriteLine(audioSample);
+ }
+
+ if (stepCount % 1000 == 0 && false)
{
- double crankDeg = cycleCrankAngle * 180.0 / Math.PI;
- double volCC = cylinder.Volume * 1e6; // cc
- Console.WriteLine($"{time,5:F3} {crankDeg,7:F1}° {volCC,5:F1} {massFlow,14:E4} {combustionCount,4} {misfireCount,4}");
+ Console.WriteLine($"{time,5:F3} {crankshaft.AngularVelocity*60/(2*Math.PI),5:F0} RPM " +
+ $"Thr:{Throttle:F2} P_target:{targetPressure/101325:F1} bar " +
+ $"mflow:{massFlow,14:E4} Comb#{engineCyl.CombustionCount} Mis#{engineCyl.MisfireCount}");
}
return audioSample;
@@ -285,7 +125,6 @@ namespace FluidSim.Core
const float T_hot = 1500f;
const float T_cold = 0f;
const float R = 287.05f;
-
float deltaHot = T_hot - T_ambient;
float deltaCold = T_ambient - T_cold;
@@ -303,12 +142,12 @@ namespace FluidSim.Core
var cylRect = new RectangleShape(new Vector2f(cylW, cylH));
cylRect.Position = new Vector2f(40f, centerY - cylH / 2f);
- double tempCyl = cylinder.Temperature; // Volume0D now has Temperature
+ double tempCyl = engineCyl.Cylinder.Temperature;
float tnCyl = NormaliseTemperature(tempCyl);
- byte redCyl = (byte)(tnCyl > 0 ? 255 * tnCyl : 0);
- byte blueCyl = (byte)(tnCyl < 0 ? -255 * tnCyl : 0);
- byte greenCyl = (byte)(255 * (1 - Math.Abs(tnCyl)));
- cylRect.FillColor = new Color(redCyl, greenCyl, blueCyl);
+ byte rC = (byte)(tnCyl > 0 ? 255 * tnCyl : 0);
+ byte bC = (byte)(tnCyl < 0 ? -255 * tnCyl : 0);
+ byte gC = (byte)(255 * (1 - Math.Abs(tnCyl)));
+ cylRect.FillColor = new Color(rC, gC, bC);
target.Draw(cylRect);
int n = exhaustPipe.GetCellCount();
@@ -317,7 +156,7 @@ namespace FluidSim.Core
float dx = pipeLen / (n - 1);
float baseRadius = 20f;
var vertices = new Vertex[n * 2];
- float ambientPressure = 101325f;
+ float ambPress = 101325f;
for (int i = 0; i < n; i++)
{
@@ -326,7 +165,7 @@ namespace FluidSim.Core
double rho = exhaustPipe.GetCellDensity(i);
double T = p / (rho * R);
- float r = baseRadius * 0.1f * (float)(1.0 + (p - ambientPressure) / ambientPressure);
+ float r = baseRadius * 0.3f * (float)(1.0 + (p - ambPress) / ambPress);
if (r < 2f) r = 2f;
float tn = NormaliseTemperature(T);