diff --git a/Components/Atmosphere.cs b/Components/Atmosphere.cs
index 1558240..8f07602 100644
--- a/Components/Atmosphere.cs
+++ b/Components/Atmosphere.cs
@@ -37,7 +37,7 @@ namespace FluidSim.Components
Port.Density = Density;
Port.Temperature = Temperature;
Port.SpecificEnthalpy = SpecificEnthalpy;
- // MassFlowRate is set by the solver connector
+ Port.AirFraction = 1.0;
}
}
}
\ No newline at end of file
diff --git a/Components/Cylinder.cs b/Components/Cylinder.cs
index 79238cb..3fc5e8c 100644
--- a/Components/Cylinder.cs
+++ b/Components/Cylinder.cs
@@ -6,7 +6,6 @@ namespace FluidSim.Components
{
public class Cylinder : IComponent
{
- // Public ports
public Port IntakePort { get; }
public Port ExhaustPort { get; }
public Crankshaft Crankshaft { get; }
@@ -20,7 +19,7 @@ namespace FluidSim.Components
public double ConRodLength { get; }
public double CompressionRatio { get; }
- // Valve timings (degrees, 0 = TDC compression, 720° full cycle)
+ // Valve timings
public double IVO { get; }
public double IVC { get; }
public double EVO { get; }
@@ -31,51 +30,48 @@ namespace FluidSim.Components
public double MaxExhaustArea { get; set; } = 0.0005;
// Ignition and combustion
- public double SparkAdvance { get; set; } = 20.0; // °BTDC
+ public double SparkAdvance { get; set; } = 20.0;
public double WiebeA { get; set; } = 5.0;
public double WiebeM { get; set; } = 2.0;
- public double WiebeDuration { get; set; } = 60.0; // degrees
- public double WiebeStart { get; set; } = 5.0; // degrees after spark
+ public double WiebeDuration { get; set; } = 60.0;
+ public double WiebeStart { get; set; } = 5.0;
// Fuel
public double StoichiometricAFR { get; set; } = 14.7;
- public double FuelLowerHeatingValue { get; set; } = 44e6; // J/kg
+ public double FuelLowerHeatingValue { get; set; } = 44e6;
// Heat loss
- public double CylinderWallArea { get; set; } = 0.02; // m²
- public double HeatTransferCoefficient { get; set; } = 100.0; // W/(m²·K)
- public double AmbientTemperature { get; set; } = 300.0; // K
+ public double CylinderWallArea { get; set; } = 0.02;
+ public double HeatTransferCoefficient { get; set; } = 100.0;
+ public double AmbientTemperature { get; set; } = 300.0;
- // State (public for drawing)
+ // State
public double Volume => cylinderVolume;
public double Pressure => (Gamma - 1.0) * cylinderEnergy / Math.Max(cylinderVolume, 1e-12);
public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
- public double Density => cylinderMass / Math.Max(cylinderVolume, 1e-12);
- public double Mass => cylinderMass;
+ public double Density => Mass / Math.Max(cylinderVolume, 1e-12);
+ public double Mass => _airMass + _exhaustMass;
+ public double AirFraction => _airMass / Math.Max(Mass, 1e-12);
public double PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
private double cylinderVolume;
- private double cylinderMass;
private double cylinderEnergy;
+ private double _airMass;
+ private double _exhaustMass;
private double trappedAirMass;
private double fuelMass;
- private double burnFraction; // 0–1
+ private double burnFraction;
private bool combustionActive;
private bool fuelInjected;
- // --- Debounce flag: allows combustion only below a certain temperature ---
- private bool _canCombust = true;
- private const double CombustionEnableTemperature = 800.0; // K – must cool below this to re‑arm
-
private const double Gamma = 1.4;
private const double GasConstant = 287.0;
- // Absolute safety limits
- private const double MaxPressurePa = 200e5; // 200 bar
- private const double MaxTemperatureK = 3500.0; // 3500 K
+ private const double MaxPressurePa = 200e5;
+ private const double MaxTemperatureK = 3500.0;
public Cylinder(double bore, double stroke, double conRodLength, double compressionRatio,
- double ivo, double ivc, double evo, double evc, double initialRPM = 1000)
+ double ivo, double ivc, double evo, double evc, Crankshaft crankshaft)
{
Bore = bore;
Stroke = stroke;
@@ -86,10 +82,12 @@ namespace FluidSim.Components
EVO = evo;
EVC = evc;
- Crankshaft = new Crankshaft(initialRPM);
+ Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
cylinderVolume = clearanceVolume;
- cylinderMass = 1.225 * clearanceVolume;
+ double initRho = 1.225;
+ _airMass = initRho * clearanceVolume;
+ _exhaustMass = 0.0;
cylinderEnergy = 101325.0 * clearanceVolume / (Gamma - 1.0);
IntakePort = new Port { Owner = this };
@@ -97,13 +95,10 @@ namespace FluidSim.Components
_ports = new[] { IntakePort, ExhaustPort };
}
- // Derived volumes
private double SweptVolume => Math.PI * 0.25 * Bore * Bore * Stroke;
private double clearanceVolume => SweptVolume / (CompressionRatio - 1.0);
private double CrankRadius => Stroke / 2.0;
private double Obliquity => CrankRadius / ConRodLength;
-
- // Crank angle in degrees (0‑720)
private double CrankDeg => (Crankshaft.CrankAngle % (4.0 * Math.PI)) * 180.0 / Math.PI % 720.0;
public double ComputeVolume(double thetaRad)
@@ -125,7 +120,6 @@ namespace FluidSim.Components
{
double deg = thetaDeg % 720.0;
if (deg < 0) deg += 720.0;
-
if (deg >= opens && deg <= closes)
{
double half = (closes - opens) * 0.5;
@@ -151,45 +145,57 @@ namespace FluidSim.Components
double crankAngleRad = Crankshaft.CrankAngle;
cylinderVolume = ComputeVolume(crankAngleRad);
- // Volume work (done BY gas, positive when expanding)
double dV = cylinderVolume - prevVolume;
+
+ // ---- Piston torque ----
+ double pRel = Pressure - 101325.0; // relative to ambient
+ double sinTh = Math.Sin(crankAngleRad);
+ double cosTh = Math.Cos(crankAngleRad);
+ double term = Math.Sqrt(1.0 - Obliquity * Obliquity * sinTh * sinTh);
+ double dxdtheta = CrankRadius * sinTh * (1.0 + Obliquity * cosTh / term);
+ double pistonArea = Math.PI * 0.25 * Bore * Bore;
+ double torque = pRel * pistonArea * dxdtheta;
+ Crankshaft.AddTorque(torque);
+
+ // Volume work (done BY gas, positive when expanding)
cylinderEnergy -= Pressure * dV;
double prevDeg = Crankshaft.PreviousAngle * 180.0 / Math.PI % 720.0;
double currDeg = crankAngleRad * 180.0 / Math.PI % 720.0;
- // ----- Intake closing: capture trapped air mass and compute fuel -----
+ // Intake closing: capture trapped air mass (air only!)
if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
{
- trappedAirMass = cylinderMass;
+ trappedAirMass = _airMass;
fuelMass = trappedAirMass / StoichiometricAFR;
fuelInjected = true;
}
- // ----- Spark ignition (once per cycle, only if canCombust) -----
+ // Spark
double sparkAngle = 0.0 - SparkAdvance;
if (sparkAngle < 0) sparkAngle += 720.0;
-
bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) ||
(prevDeg > sparkAngle + 360.0 && currDeg < sparkAngle);
- if (crossedSpark && !combustionActive && fuelInjected && _canCombust)
+ if (crossedSpark && !combustionActive && fuelInjected)
{
combustionActive = true;
burnFraction = 0.0;
}
- // ----- Combustion progress -----
+ // Combustion progress
if (combustionActive)
{
double angleSinceSpark = currDeg - sparkAngle;
if (angleSinceSpark < 0) angleSinceSpark += 720.0;
double newFraction = Wiebe(angleSinceSpark);
-
if (newFraction >= 1.0 || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance))
{
newFraction = 1.0;
combustionActive = false;
- _canCombust = false; // require cool‑down before next ignition
+ // All gas becomes exhaust
+ double totalMass = _airMass + _exhaustMass;
+ _airMass = 0.0;
+ _exhaustMass = totalMass;
}
double dFraction = newFraction - burnFraction;
@@ -197,18 +203,12 @@ namespace FluidSim.Components
{
double dQ = fuelMass * FuelLowerHeatingValue * dFraction;
cylinderEnergy += dQ;
- cylinderMass += fuelMass * dFraction;
+ _exhaustMass += fuelMass * dFraction; // burning fuel adds to exhaust
burnFraction = newFraction;
}
}
- // ----- Re‑arm combustion if temperature has dropped low enough -----
- if (!combustionActive && !_canCombust && Temperature < CombustionEnableTemperature)
- {
- _canCombust = true;
- }
-
- // ----- Heat loss to cylinder walls -----
+ // Heat loss
double dQ_loss = HeatTransferCoefficient * CylinderWallArea *
(Temperature - AmbientTemperature) * dt;
cylinderEnergy -= dQ_loss;
@@ -216,39 +216,46 @@ namespace FluidSim.Components
// Update port states
double p = Pressure, rho = Density, T = Temperature;
double h = Gamma / (Gamma - 1.0) * p / Math.Max(rho, 1e-12);
+ double af = AirFraction;
+
IntakePort.Pressure = p;
IntakePort.Density = rho;
IntakePort.Temperature = T;
IntakePort.SpecificEnthalpy = h;
+ IntakePort.AirFraction = af;
ExhaustPort.Pressure = p;
ExhaustPort.Density = rho;
ExhaustPort.Temperature = T;
ExhaustPort.SpecificEnthalpy = h;
+ ExhaustPort.AirFraction = af;
}
public void UpdateState(double dt)
{
- double dm = 0.0;
- double dE = 0.0;
+ double dmAir = 0.0, dmExhaust = 0.0, dE = 0.0;
foreach (var port in _ports)
{
- dm += port.MassFlowRate * dt;
- dE += port.MassFlowRate * port.SpecificEnthalpy * dt;
+ double mdot = port.MassFlowRate;
+ double af = mdot >= 0 ? port.AirFraction : AirFraction;
+ dmAir += mdot * af * dt;
+ dmExhaust += mdot * (1.0 - af) * dt;
+ dE += mdot * port.SpecificEnthalpy * dt;
}
- cylinderMass += dm;
+ _airMass += dmAir;
+ _exhaustMass += dmExhaust;
cylinderEnergy += dE;
double V = Math.Max(cylinderVolume, 1e-12);
- // --- Absolute pressure & temperature clamps ---
+ // Safety clamps
double currentP = (Gamma - 1.0) * cylinderEnergy / V;
if (currentP > MaxPressurePa)
cylinderEnergy = MaxPressurePa * V / (Gamma - 1.0);
- double currentRho = cylinderMass / V;
+ double currentRho = (_airMass + _exhaustMass) / V;
double currentT = currentP / Math.Max(currentRho * GasConstant, 1e-12);
if (currentT > MaxTemperatureK)
{
@@ -256,10 +263,11 @@ namespace FluidSim.Components
cylinderEnergy = pAtTlimit * V / (Gamma - 1.0);
}
- // Existing safeguards
- if (cylinderMass < 1e-9)
+ double totalMass = _airMass + _exhaustMass;
+ if (totalMass < 1e-9)
{
- cylinderMass = 1e-9;
+ _airMass = 1e-9;
+ _exhaustMass = 0.0;
cylinderEnergy = 101325.0 * V / (Gamma - 1.0);
}
else if (cylinderEnergy < 0.0)
@@ -267,8 +275,8 @@ namespace FluidSim.Components
cylinderEnergy = 101325.0 * V / (Gamma - 1.0);
}
- if (cylinderMass < 0.0) cylinderMass = 1e-9;
- if (cylinderEnergy < 0.0) cylinderEnergy = 101325.0 * V / (Gamma - 1.0);
+ if (_airMass < 0.0) _airMass = 0.0;
+ if (_exhaustMass < 0.0) _exhaustMass = 0.0;
}
}
}
\ No newline at end of file
diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
index 4b0a122..a0bbcd3 100644
--- a/Components/Pipe1D.cs
+++ b/Components/Pipe1D.cs
@@ -7,10 +7,10 @@ namespace FluidSim.Components
///
/// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme.
/// Ghost states are set externally via SetGhostLeft/Right; they are always required.
+ /// Now includes a passive scalar for air mass fraction.
///
public class Pipe1D : IComponent
{
- // ---------- Compile‑time profiling flag ----------
public const bool EnableDetailedProfiling = false; // set to false in release builds
public Port PortA { get; }
@@ -36,10 +36,11 @@ namespace FluidSim.Components
private readonly double _gamma = 1.4;
private double[] _rho, _rhou, _E;
- private double[] _fluxM, _fluxP, _fluxE; // flux at cell faces (0.._n) – kept for possible external use, not used internally anymore
+ private double[] _Y; // air mass fraction
+ private double[] _fluxM, _fluxP, _fluxE;
- private double _rhoGhostL, _uGhostL, _pGhostL;
- private double _rhoGhostR, _uGhostR, _pGhostR;
+ private double _rhoGhostL, _uGhostL, _pGhostL, _YGhostL;
+ private double _rhoGhostR, _uGhostR, _pGhostR, _YGhostR;
private bool _ghostLValid, _ghostRValid;
private double _laminarCoeff;
@@ -65,7 +66,7 @@ namespace FluidSim.Components
_rho = new double[_n];
_rhou = new double[_n];
_E = new double[_n];
-
+ _Y = new double[_n];
_fluxM = new double[_n + 1];
_fluxP = new double[_n + 1];
_fluxE = new double[_n + 1];
@@ -87,16 +88,19 @@ namespace FluidSim.Components
public void UpdateState(double dt) { }
// ---------- Ghost interface ----------
- public void SetGhostLeft(double rho, double u, double p)
+ public void SetGhostLeft(double rho, double u, double p, double airFraction)
{
- _rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _ghostLValid = true;
+ _rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _YGhostL = airFraction; _ghostLValid = true;
}
- public void SetGhostRight(double rho, double u, double p)
+ public void SetGhostRight(double rho, double u, double p, double airFraction)
{
- _rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _ghostRValid = true;
+ _rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _YGhostR = airFraction; _ghostRValid = true;
}
public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; }
+ public double GetInteriorAirFractionLeft() => _Y[0];
+ public double GetInteriorAirFractionRight() => _Y[_n - 1];
+
public (double rho, double u, double p) GetInteriorStateLeft()
{
double rho = Math.Max(_rho[0], 1e-12);
@@ -172,6 +176,34 @@ namespace FluidSim.Components
t0 = t1;
}
+ // ---------- Local flux functions ----------
+ void LaxFlux(double rL, double uL, double pL, double cL,
+ double rR, double uR, double pR, double cR,
+ out double fm, out double fp, out double fe)
+ {
+ double EL = pL / (gm1 * rL) + 0.5 * uL * uL;
+ double ER = pR / (gm1 * rR) + 0.5 * uR * uR;
+ double Fm_L = rL * uL;
+ double Fp_L = rL * uL * uL + pL;
+ double Fe_L = (rL * EL + pL) * uL;
+ double Fm_R = rR * uR;
+ double Fp_R = rR * uR * uR + pR;
+ double Fe_R = (rR * ER + pR) * uR;
+ double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
+ fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rR - rL);
+ fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rR * uR - rL * uL);
+ fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (rR * ER - rL * EL);
+ }
+
+ void ScalarFlux(double rL, double uL, double cL, double YL,
+ double rR, double uR, double cR, double YR,
+ double alpha, out double fy)
+ {
+ double Fm_L = rL * uL;
+ double Fm_R = rR * uR;
+ fy = 0.5 * (Fm_L * YL + Fm_R * YR) - 0.5 * alpha * (rR * YR - rL * YL);
+ }
+
// ---------- Phase 2: Left face flux (ghostL – cell 0) ----------
double rL_ghost = Math.Max(_rhoGhostL, 1e-12);
double pL_ghost = _pGhostL;
@@ -182,6 +214,12 @@ namespace FluidSim.Components
_rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), p[0], c[0],
out double fluxM_left, out double fluxP_left, out double fluxE_left);
+ double alphaLeft = Math.Max(Math.Abs(uL_ghost) + cL_ghost,
+ Math.Abs(_rhou[0] / Math.Max(_rho[0], 1e-12)) + c[0]);
+ ScalarFlux(rL_ghost, uL_ghost, cL_ghost, _YGhostL,
+ _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), c[0], _Y[0],
+ alphaLeft, out double fluxY_left);
+
if (EnableDetailedProfiling)
{
t1 = Stopwatch.GetTimestamp();
@@ -193,6 +231,7 @@ namespace FluidSim.Components
double fluxM_prev = fluxM_left;
double fluxP_prev = fluxP_left;
double fluxE_prev = fluxE_left;
+ double fluxY_prev = fluxY_left;
for (int i = 0; i < n - 1; i++)
{
@@ -203,23 +242,31 @@ namespace FluidSim.Components
double uL = _rhou[iL] / rL;
double pL = p[iL];
double cL = c[iL];
+ double YL = _Y[iL];
double rR = Math.Max(_rho[iR], 1e-12);
double uR = _rhou[iR] / rR;
double pR = p[iR];
double cR = c[iR];
+ double YR = _Y[iR];
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR,
out double fluxM_right, out double fluxP_right, out double fluxE_right);
+ double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
+ ScalarFlux(rL, uL, cL, YL, rR, uR, cR, YR, alpha, out double fluxY_right);
+
// Update cell i
double r = _rho[i];
double ru = _rhou[i];
double E = _E[i];
+ double Y = _Y[i];
double newR = r - dt_dx * (fluxM_right - fluxM_prev);
double newRu = ru - dt_dx * (fluxP_right - fluxP_prev);
double newE = E - dt_dx * (fluxE_right - fluxE_prev);
+ double oldRhoY = r * Y;
+ double newRhoY = oldRhoY - dt_dx * (fluxY_right - fluxY_prev);
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
newRu *= dampingFactor;
@@ -234,10 +281,12 @@ namespace FluidSim.Components
_rho[i] = newR;
_rhou[i] = newRu;
_E[i] = newE;
+ _Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0);
fluxM_prev = fluxM_right;
fluxP_prev = fluxP_right;
fluxE_prev = fluxE_right;
+ fluxY_prev = fluxY_right;
}
if (EnableDetailedProfiling)
@@ -253,20 +302,31 @@ namespace FluidSim.Components
double uR_ghost = _uGhostR;
double cR_ghost = Math.Sqrt(gamma * pR_ghost / rR_ghost);
- LaxFlux(_rho[n - 1], _rhou[n - 1] / Math.Max(_rho[n - 1], 1e-12), p[n - 1], c[n - 1],
+ double rInt = _rho[n - 1];
+ double uInt = _rhou[n - 1] / Math.Max(rInt, 1e-12);
+
+ LaxFlux(rInt, uInt, p[n - 1], c[n - 1],
rR_ghost, uR_ghost, pR_ghost, cR_ghost,
out double fluxM_right_final, out double fluxP_right_final, out double fluxE_right_final);
- // Update last cell (identical to interior, but with final fluxes)
+ double alphaRight = Math.Max(Math.Abs(uInt) + c[n - 1], Math.Abs(uR_ghost) + cR_ghost);
+ ScalarFlux(rInt, uInt, c[n - 1], _Y[n - 1],
+ rR_ghost, uR_ghost, cR_ghost, _YGhostR,
+ alphaRight, out double fluxY_right_final);
+
+ // Update last cell
{
int i = n - 1;
double r = _rho[i];
double ru = _rhou[i];
double E = _E[i];
+ double Y = _Y[i];
double newR = r - dt_dx * (fluxM_right_final - fluxM_prev);
double newRu = ru - dt_dx * (fluxP_right_final - fluxP_prev);
double newE = E - dt_dx * (fluxE_right_final - fluxE_prev);
+ double oldRhoY = r * Y;
+ double newRhoY = oldRhoY - dt_dx * (fluxY_right_final - fluxY_prev);
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
newRu *= dampingFactor;
@@ -281,6 +341,7 @@ namespace FluidSim.Components
_rho[i] = newR;
_rhou[i] = newRu;
_E[i] = newE;
+ _Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0);
}
if (EnableDetailedProfiling)
@@ -295,11 +356,13 @@ namespace FluidSim.Components
PortA.Pressure = pA; PortA.Density = rhoA;
PortA.Temperature = pA / (rhoA * 287.0);
PortA.SpecificEnthalpy = gm1 / (gamma - 1.0) * pA / rhoA;
+ PortA.AirFraction = _Y[0];
(double rhoB, double uB, double pB) = GetInteriorStateRight();
PortB.Pressure = pB; PortB.Density = rhoB;
PortB.Temperature = pB / (rhoB * 287.0);
PortB.SpecificEnthalpy = gm1 / (gamma - 1.0) * pB / rhoB;
+ PortB.AirFraction = _Y[_n - 1];
if (EnableDetailedProfiling)
{
@@ -308,48 +371,6 @@ namespace FluidSim.Components
}
}
- // ---------- Local Lax‑Friedrichs flux function ----------
- private void LaxFlux(double rL, double uL, double pL, double cL,
- double rR, double uR, double pR, double cR,
- out double fm, out double fp, out double fe)
- {
- double gm1 = _gamma - 1.0;
- double EL = pL / (gm1 * rL) + 0.5 * uL * uL;
- double ER = pR / (gm1 * rR) + 0.5 * uR * uR;
- double Fm_L = rL * uL;
- double Fp_L = rL * uL * uL + pL;
- double Fe_L = (rL * EL + pL) * uL;
- double Fm_R = rR * uR;
- double Fp_R = rR * uR * uR + pR;
- double Fe_R = (rR * ER + pR) * uR;
- double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
- fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rR - rL);
- fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rR * uR - rL * uL);
- fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (rR * ER - rL * EL);
- }
-
- // Original LaxFriedrichsFlux (kept for compatibility, can be removed if unused)
- private void LaxFriedrichsFlux(double rL, double uL, double pL, double eL,
- double rR, double uR, double pR, double eR,
- out double fm, out double fp, out double fe)
- {
- double rhoL = rL, rhoR = rR;
- double EL = rhoL * eL;
- double ER = rhoR * eR;
- double Fm_L = rhoL * uL;
- double Fp_L = rhoL * uL * uL + pL;
- double Fe_L = (EL + pL) * uL;
- double Fm_R = rhoR * uR;
- double Fp_R = rhoR * uR * uR + pR;
- double Fe_R = (ER + pR) * uR;
- double cL = Math.Sqrt(_gamma * pL / rL);
- double cR = Math.Sqrt(_gamma * pR / rR);
- double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
- fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rhoR - rhoL);
- fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rhoR * uR - rhoL * uL);
- fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (ER - EL);
- }
-
private double PressureScalar(int i)
{
double rho = Math.Max(_rho[i], 1e-12);
@@ -365,6 +386,7 @@ namespace FluidSim.Components
_rho[i] = rho;
_rhou[i] = rho * u;
_E[i] = E;
+ _Y[i] = 1.0; // initially pure air
}
}
@@ -376,6 +398,7 @@ namespace FluidSim.Components
_rho[i] = rho;
_rhou[i] = rho * u;
_E[i] = E;
+ _Y[i] = 1.0;
}
public void SetCellPressure(int i, double p)
diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs
index 14103e0..35b517f 100644
--- a/Components/Volume0D.cs
+++ b/Components/Volume0D.cs
@@ -1,4 +1,4 @@
-using System;
+using System;
using System.Collections.Generic;
using FluidSim.Interfaces;
@@ -8,8 +8,9 @@ namespace FluidSim.Components
{
public List Ports { get; } = new List();
- public double Mass { get; set; } // made public setter
- public double InternalEnergy { get; set; } // made public setter
+ private double _airMass;
+ private double _exhaustMass;
+ public double InternalEnergy { get; set; }
public double Volume { get; set; }
public double Dvdt { get; set; }
public double Gamma { get; set; } = 1.4;
@@ -18,6 +19,8 @@ namespace FluidSim.Components
public double AmbientPressure { get; set; } = 101325.0;
// Derived quantities
+ public double Mass => _airMass + _exhaustMass;
+ public double AirFraction => _airMass / Math.Max(Mass, 1e-12);
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);
@@ -32,7 +35,8 @@ namespace FluidSim.Components
Dvdt = 0.0;
double rho0 = initialPressure / (GasConstant * initialTemperature);
- Mass = rho0 * Volume;
+ _airMass = rho0 * Volume; // starts with all air
+ _exhaustMass = 0.0;
InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
}
@@ -43,44 +47,53 @@ namespace FluidSim.Components
port.Density = Density;
port.Temperature = Temperature;
port.SpecificEnthalpy = SpecificEnthalpy;
+ port.AirFraction = AirFraction;
Ports.Add(port);
return port;
}
- ///
- /// Set the pressure to a specific value while keeping the current temperature constant.
- /// Updates Mass and InternalEnergy accordingly.
- ///
public void SetPressure(double pressure, double? temperature = null)
{
double V = Math.Max(Volume, 1e-12);
double T = temperature ?? Temperature;
double rho = pressure / (GasConstant * T);
- Mass = rho * V;
+ double totalMass = rho * V;
+ // Keep current air fraction when setting pressure?
+ double af = AirFraction;
+ _airMass = totalMass * af;
+ _exhaustMass = totalMass * (1.0 - af);
InternalEnergy = pressure * V / (Gamma - 1.0);
}
public void UpdateState(double dt)
{
- double totalMdot = 0.0;
+ double totalMdotAir = 0.0;
+ double totalMdotExhaust = 0.0;
double totalEdot = 0.0;
foreach (var port in Ports)
{
- totalMdot += port.MassFlowRate;
- totalEdot += port.MassFlowRate * port.SpecificEnthalpy;
+ double mdot = port.MassFlowRate; // positive INTO volume
+ double af = mdot >= 0 ? port.AirFraction : AirFraction; // inflow: use port's fraction; outflow: well-mixed
+ totalMdotAir += mdot * af;
+ totalMdotExhaust += mdot * (1.0 - af);
+ totalEdot += mdot * port.SpecificEnthalpy;
}
- double dm = totalMdot * dt;
+ double dAir = totalMdotAir * dt;
+ double dExhaust = totalMdotExhaust * dt;
double dE = totalEdot * dt - Pressure * Dvdt * dt;
- Mass += dm;
+ _airMass += dAir;
+ _exhaustMass += dExhaust;
InternalEnergy += dE;
double V = Math.Max(Volume, 1e-12);
- if (Mass < 1e-9)
+ double totalMass = _airMass + _exhaustMass;
+ if (totalMass < 1e-9)
{
- Mass = 1e-9;
+ _airMass = 1e-9;
+ _exhaustMass = 0.0;
InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
}
else if (InternalEnergy < 0.0)
@@ -88,16 +101,17 @@ namespace FluidSim.Components
InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
}
- if (Mass < 0.0) Mass = 1e-9;
- if (InternalEnergy < 0.0) InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
+ if (_airMass < 0.0) _airMass = 0.0;
+ if (_exhaustMass < 0.0) _exhaustMass = 0.0;
- double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy;
+ double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afrac = AirFraction;
foreach (var port in Ports)
{
port.Pressure = p;
port.Density = rho;
port.Temperature = T;
port.SpecificEnthalpy = h;
+ port.AirFraction = afrac;
}
}
diff --git a/Core/Junction.cs b/Core/Junction.cs
deleted file mode 100644
index 3972b90..0000000
--- a/Core/Junction.cs
+++ /dev/null
@@ -1,228 +0,0 @@
-using System;
-using System.Collections.Generic;
-using FluidSim.Components;
-
-namespace FluidSim.Core
-{
- ///
- /// Zero‑dimensional junction connecting multiple pipe ends.
- /// The coupling conditions are mass conservation and equality of
- /// stagnation enthalpy (Bernoulli invariant) for all branches,
- /// following Reigstad (2014, 2015). A root‑finding method (Brent)
- /// solves for the common junction pressure.
- ///
- public class Junction
- {
- public struct Branch
- {
- public Pipe1D Pipe;
- public bool IsLeftEnd;
- }
-
- private readonly List _branches = new List();
- public IReadOnlyList Branches => _branches;
-
- // Last resolved state (for audio / monitoring)
- public double LastJunctionPressure { get; private set; }
- public double[] LastBranchMassFlows { get; private set; } = Array.Empty();
-
- public Junction() { }
-
- public void AddBranch(Pipe1D pipe, bool isLeftEnd)
- {
- _branches.Add(new Branch { Pipe = pipe, IsLeftEnd = isLeftEnd });
- }
-
- ///
- /// Solve the junction for one sub‑step. Uses Brent's method to find
- /// the pressure p* that satisfies sum(mdot) = 0 with stagnation enthalpy equality.
- ///
- public void Resolve(double dtSub)
- {
- int nb = _branches.Count;
- if (nb < 2)
- throw new InvalidOperationException("Junction requires at least 2 branches.");
-
- // Gather interior states and areas
- var rho = new double[nb];
- var u = new double[nb];
- var p = new double[nb];
- var area = new double[nb];
- var isLeft = new bool[nb];
- double gamma = 1.4;
-
- double pMin = double.MaxValue, pMax = double.MinValue;
- for (int i = 0; i < nb; i++)
- {
- var branch = _branches[i];
- (double ri, double ui, double pi) = branch.IsLeftEnd
- ? branch.Pipe.GetInteriorStateLeft()
- : branch.Pipe.GetInteriorStateRight();
- rho[i] = ri; u[i] = ui; p[i] = pi;
- area[i] = branch.Pipe.Area;
- isLeft[i] = branch.IsLeftEnd;
-
- if (pi < pMin) pMin = pi;
- if (pi > pMax) pMax = pi;
- }
-
- // We solve for pStar that makes totalMassFlow(pStar) = 0.
- // The function: totalMassFlow = sum( sign_i * rhoStar_i * uStar_i * A_i )
- // where for each branch:
- // - Riemann invariant: J = u + 2c/(γ-1) for right end, J = u - 2c/(γ-1) for left end.
- // - uStar = J ∓ 2cStar/(γ-1) (depending on direction)
- // - Isentropic relation: rhoStar = rho_i * (pStar / p_i)^{1/γ}
- // - cStar = sqrt(γ pStar / rhoStar)
- // We require stagnation enthalpy equality: h0 = h + u^2/2 = constant across junction.
- // Hence for each branch we compute the specific total enthalpy:
- // hStar = (γ/(γ-1)) * pStar/rhoStar, h0_star = hStar + 0.5 uStar^2.
- // We enforce that all h0_star are equal. Mass conservation then determines pStar.
- // This is a scalar root‑finding problem.
-
- // Bracket the solution: pressure must lie between min and max of branch pressures (expanded a bit)
- double a = Math.Max(100.0, pMin * 0.1);
- double b = Math.Min(1e7, pMax * 10.0);
- if (a >= b) { a = 100.0; b = 1e7; }
-
- Func f = pStar =>
- {
- double totalMdot = 0.0;
- double h0Ref = 0.0;
- bool first = true;
- for (int i = 0; i < nb; i++)
- {
- double g = gamma;
- double gm1 = g - 1.0;
- double rhoI = rho[i], uI = u[i], pI = p[i];
- double cI = Math.Sqrt(g * pI / rhoI);
- double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1;
-
- double pratio = Math.Max(pStar / pI, 1e-6);
- double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g);
- double cStar = Math.Sqrt(g * pStar / rhoStar);
- double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1;
-
- double hStar = (g / gm1) * pStar / rhoStar;
- double h0 = hStar + 0.5 * uStar * uStar;
-
- if (first)
- {
- h0Ref = h0;
- first = false;
- }
- else
- {
- // Equality of stagnation enthalpy: ideally h0 == h0Ref.
- // We incorporate a penalty to enforce this.
- }
-
- // Mass flow into junction: sign convention = positive if fluid leaves pipe into junction.
- double sign = isLeft[i] ? -1.0 : 1.0; // left end: positive u is into pipe, so into junction is -u
- double mdot_i = sign * rhoStar * uStar * area[i];
- totalMdot += mdot_i;
- }
-
- // Additional term to enforce equal stagnation enthalpies? For simplicity, we only enforce mass conservation here,
- // because with the Riemann invariants and a common pressure, the stagnation enthalpies are automatically equal
- // if the junction is isentropic? Actually, with a common pressure and isentropic relations from each branch,
- // each branch has its own entropy (p/ρ^γ = const), so h0 may differ. The correct condition is mass conservation + equality of h0.
- // To solve both, we would need to vary pStar and a common h0? In Reigstad's formulation, the system yields
- // mass conservation as the determinant, and pStar is found from that equation, with the assumption that the junction
- // itself does not introduce entropy. The typical implementation uses the Riemann invariants and mass conservation only.
- // We'll stick to mass conservation for now.
- return totalMdot;
- };
-
- double pStar = BrentsMethod(f, a, b, 1e-6, 100);
- LastJunctionPressure = pStar;
- LastBranchMassFlows = new double[nb];
-
- // Apply ghost states and record mass flows
- for (int i = 0; i < nb; i++)
- {
- double g = gamma, gm1 = g - 1.0;
- double rhoI = rho[i], uI = u[i], pI = p[i];
- double cI = Math.Sqrt(g * pI / rhoI);
- double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1;
-
- double pratio = Math.Max(pStar / pI, 1e-6);
- double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g);
- double cStar = Math.Sqrt(g * pStar / rhoStar);
- double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1;
-
- double sign = isLeft[i] ? -1.0 : 1.0;
- double mdot = sign * rhoStar * uStar * area[i];
- LastBranchMassFlows[i] = mdot;
-
- if (isLeft[i])
- _branches[i].Pipe.SetGhostLeft(rhoStar, uStar, pStar);
- else
- _branches[i].Pipe.SetGhostRight(rhoStar, uStar, pStar);
- }
- }
-
- /// Simple Brent's method root finder.
- private static double BrentsMethod(Func f, double a, double b, double tol, int maxIter)
- {
- double fa = f(a), fb = f(b);
- if (fa * fb >= 0)
- return (a + b) / 2.0; // fallback
-
- double c = a, fc = fa;
- double d = b - a, e = d;
-
- for (int iter = 0; iter < maxIter; iter++)
- {
- if (Math.Abs(fc) < Math.Abs(fb))
- {
- a = b; b = c; c = a;
- fa = fb; fb = fc; fc = fa;
- }
- double tol1 = 2 * double.Epsilon * Math.Abs(b) + 0.5 * tol;
- double xm = 0.5 * (c - b);
- if (Math.Abs(xm) <= tol1 || fb == 0.0)
- return b;
-
- if (Math.Abs(e) >= tol1 && Math.Abs(fa) > Math.Abs(fb))
- {
- double s = fb / fa;
- double p, q;
- if (a == c)
- {
- p = 2.0 * xm * s;
- q = 1.0 - s;
- }
- else
- {
- q = fa / fc;
- double r = fb / fc;
- p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
- q = (q - 1.0) * (r - 1.0) * (s - 1.0);
- }
- if (p > 0) q = -q; else p = -p;
- s = e; e = d;
- if (2.0 * p < 3.0 * xm * q - Math.Abs(tol1 * q) && p < Math.Abs(0.5 * s * q))
- {
- d = p / q;
- }
- else
- {
- d = xm; e = d;
- }
- }
- else
- {
- d = xm; e = d;
- }
-
- a = b; fa = fb;
- if (Math.Abs(d) > tol1)
- b += d;
- else
- b += Math.Sign(xm) * tol1;
- fb = f(b);
- }
- return b;
- }
- }
-}
\ No newline at end of file
diff --git a/Core/OpenEndLink.cs b/Core/OpenEndLink.cs
index ca40d4a..ba57d13 100644
--- a/Core/OpenEndLink.cs
+++ b/Core/OpenEndLink.cs
@@ -9,6 +9,7 @@ namespace FluidSim.Core
/// from the isentropic expansion to ambient pressure, using the pipe's entropy,
/// and the outgoing Riemann invariant. This avoids a density jump at flow reversal.
/// Supersonic outflow extrapolates the interior state.
+ /// Now includes air fraction tracking: incoming air is fresh (AF=1), outgoing uses interior pipe AF.
///
public class OpenEndLink
{
@@ -34,34 +35,34 @@ namespace FluidSim.Core
? Pipe.GetInteriorStateLeft()
: Pipe.GetInteriorStateRight();
+ double airFracInt = IsLeftEnd
+ ? Pipe.GetInteriorAirFractionLeft()
+ : Pipe.GetInteriorAirFractionRight();
+
double gamma = Gamma;
double gm1 = gamma - 1.0;
double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12));
double pAmb = AmbientPressure;
// Riemann invariants
- double J_plus = uInt + 2.0 * cInt / gm1;
+ double J_plus = uInt + 2.0 * cInt / gm1;
double J_minus = uInt - 2.0 * cInt / gm1;
- double rhoGhost, uGhost, pGhost;
+ double rhoGhost, uGhost, pGhost, airFracGhost;
// ---- Subsonic branch (used for both outflow and inflow) ----
- // Isentropic expansion to ambient pressure using pipe's entropy
double s = pInt / Math.Pow(rhoInt, gamma); // entropy constant
double rhoIso = Math.Pow(pAmb / s, 1.0 / gamma);
double cIso = Math.Sqrt(gamma * pAmb / Math.Max(rhoIso, 1e-12));
double uIso = IsLeftEnd
? (J_minus + 2.0 * cIso / gm1)
- : (J_plus - 2.0 * cIso / gm1);
+ : (J_plus - 2.0 * cIso / gm1);
- // Check for supersonic outflow: if the isentropic velocity exceeds the speed of sound,
- // the flow is supersonic and we extrapolate the interior state.
+ // Check for supersonic outflow
bool supersonic = IsLeftEnd
- ? (uInt <= -cInt) // left end: outflow is when u < -c
- : (uInt >= cInt); // right end: outflow is when u > c
+ ? (uInt <= -cInt)
+ : (uInt >= cInt);
- // Extra check: if the isentropic velocity is supersonic in the outflow direction,
- // also treat as supersonic (this can happen when the interior pressure is very high).
if (!supersonic)
{
if (IsLeftEnd)
@@ -74,27 +75,42 @@ namespace FluidSim.Core
{
// Supersonic outflow – extrapolate interior
rhoGhost = rhoInt;
- uGhost = uInt;
- pGhost = pInt;
+ uGhost = uInt;
+ pGhost = pInt;
+ airFracGhost = airFracInt; // whatever is leaving
}
else
{
- // Subsonic flow – use the isentropic state
+ // Subsonic flow – use isentropic state to ambient
rhoGhost = rhoIso;
- uGhost = uIso;
- pGhost = pAmb;
+ uGhost = uIso;
+ pGhost = pAmb;
+
+ // Determine if inflow or outflow
+ bool isInflow = IsLeftEnd ? (uIso >= 0) : (uIso <= 0); // positive u means into pipe from left end? Wait: left end: u>0 means flow to the right, into pipe. Right end: u>0 means flow to the right, out of pipe. Let's use mass flow sign later.
+ // More straightforward: if using the isentropic state, the ghost velocity direction indicates flow. For inflow (ambient to pipe), airFraction = 1.0; for outflow, airFraction = interior's AF.
+ if ((IsLeftEnd && uIso >= 0) || (!IsLeftEnd && uIso <= 0))
+ {
+ // Inflow (ambient enters pipe)
+ airFracGhost = 1.0;
+ }
+ else
+ {
+ // Outflow (pipe exits to ambient)
+ airFracGhost = airFracInt;
+ }
}
// Apply ghost to pipe
if (IsLeftEnd)
- Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost);
+ Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost, airFracGhost);
else
- Pipe.SetGhostRight(rhoGhost, uGhost, pGhost);
+ Pipe.SetGhostRight(rhoGhost, uGhost, pGhost, airFracGhost);
// Mass flow out of the pipe (positive = leaving)
double area = Pipe.Area;
double mdot = rhoGhost * uGhost * area;
- if (IsLeftEnd) mdot = -mdot; // left end: positive u is into pipe, so out is -u
+ if (IsLeftEnd) mdot = -mdot; // left end: positive u is into pipe, outward flow is -u
LastMassFlowRate = mdot;
LastFaceDensity = rhoGhost;
LastFaceVelocity = uGhost;
diff --git a/Core/OrificeLink.cs b/Core/OrificeLink.cs
index bc4cc4d..a416e6f 100644
--- a/Core/OrificeLink.cs
+++ b/Core/OrificeLink.cs
@@ -46,16 +46,21 @@ namespace FluidSim.Core
double volRho = VolumePort.Density;
double volT = VolumePort.Temperature;
double volH = VolumePort.SpecificEnthalpy;
+ double volAF = VolumePort.AirFraction;
(double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd
? Pipe.GetInteriorStateLeft()
: Pipe.GetInteriorStateRight();
double pipeT = pipeP / Math.Max(pipeRho * 287.0, 1e-12);
+ double pipeAF = IsPipeLeftEnd
+ ? Pipe.GetInteriorAirFractionLeft()
+ : Pipe.GetInteriorAirFractionRight();
+
double gamma = 1.4;
double R = 287.0;
- // ---- Steady‑state nozzle solution (gives correct exit state) ----
+ // ---- Steady‑state nozzle solution ----
double mdotSS; // positive = volume → pipe
double rhoFace0, uFace0, pFace0;
if (volP >= pipeP)
@@ -71,15 +76,6 @@ namespace FluidSim.Core
mdotSS = -mdotUpToDown;
}
- // ====== Hard physical cap: max sonic flow × 1.1 ======
- double upRho = mdotSS >= 0 ? volRho : pipeRho;
- double upT = mdotSS >= 0 ? volT : pipeT;
- double upC = Math.Sqrt(gamma * R * upT);
- double maxFlow = upRho * upC * area * 1.1;
- if (Math.Abs(mdotSS) > maxFlow)
- mdotSS = Math.Sign(mdotSS) * maxFlow;
- // ====================================================
-
// ---- Dynamic update ----
if (UseInertance)
{
@@ -102,21 +98,38 @@ namespace FluidSim.Core
if (_mdot > maxOut) _mdot = maxOut;
}
- // ---- Ghost state ----
+ // ---- Ghost state with air fraction ----
double rhoFace = _mdot >= 0 ? volRho : pipeRho;
double pFace = pFace0;
double mdotMag = Math.Abs(_mdot);
double uFace = mdotMag / (rhoFace * area);
+ // Determine air fraction for ghost and for volume port
+ double airFracGhost; // air fraction of ghost cell (at pipe end)
+ double airFracForVolume; // if flow reverses into volume, this is the air fraction entering volume
+
+ if (_mdot >= 0) // volume → pipe
+ {
+ airFracGhost = volAF;
+ // Flow enters pipe; no need to set volume's air fraction (port already has its own)
+ airFracForVolume = volAF; // unused
+ }
+ else // pipe → volume
+ {
+ airFracGhost = pipeAF;
+ airFracForVolume = pipeAF;
+ VolumePort.AirFraction = airFracForVolume;
+ }
+
if (IsPipeLeftEnd)
uFace = _mdot >= 0 ? uFace : -uFace;
else
uFace = _mdot >= 0 ? -uFace : uFace;
if (IsPipeLeftEnd)
- Pipe.SetGhostLeft(rhoFace, uFace, pFace);
+ Pipe.SetGhostLeft(rhoFace, uFace, pFace, airFracGhost);
else
- Pipe.SetGhostRight(rhoFace, uFace, pFace);
+ Pipe.SetGhostRight(rhoFace, uFace, pFace, airFracGhost);
// Store results (positive = into volume)
LastMassFlowRate = -_mdot;
@@ -145,9 +158,9 @@ namespace FluidSim.Core
: Pipe.GetInteriorStateRight();
if (IsPipeLeftEnd)
- Pipe.SetGhostLeft(rInt, -uInt, pInt);
+ Pipe.SetGhostLeft(rInt, -uInt, pInt, IsPipeLeftEnd ? Pipe.GetInteriorAirFractionLeft() : Pipe.GetInteriorAirFractionRight());
else
- Pipe.SetGhostRight(rInt, -uInt, pInt);
+ Pipe.SetGhostRight(rInt, -uInt, pInt, IsPipeLeftEnd ? Pipe.GetInteriorAirFractionLeft() : Pipe.GetInteriorAirFractionRight());
LastMassFlowRate = 0.0;
LastFaceDensity = rInt;
diff --git a/Core/OutdoorExhaustReverb.cs b/Core/OutdoorExhaustReverb.cs
index 23dc997..517e4fb 100644
--- a/Core/OutdoorExhaustReverb.cs
+++ b/Core/OutdoorExhaustReverb.cs
@@ -19,7 +19,7 @@ namespace FluidSim.Core
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 Feedback { get; set; } = 0.55f; // safe range 0.7‑0.9
public float DampingFreq { get; set; } = 6000f; // Hz
public OutdoorExhaustReverb(int sampleRate)
diff --git a/Core/Solver.cs b/Core/Solver.cs
index 5ef78a6..ab2e8f4 100644
--- a/Core/Solver.cs
+++ b/Core/Solver.cs
@@ -11,7 +11,6 @@ namespace FluidSim.Core
{
private readonly List _components = new();
private readonly List _orificeLinks = new();
- private readonly List _junctions = new();
private readonly List _openEndLinks = new();
private double _dt;
@@ -37,7 +36,6 @@ namespace FluidSim.Core
public void AddComponent(IComponent component) => _components.Add(component);
public void AddOrificeLink(OrificeLink link) => _orificeLinks.Add(link);
- public void AddJunction(Junction junction) => _junctions.Add(junction);
public void AddOpenEndLink(OpenEndLink link) => _openEndLinks.Add(link);
public void Step()
@@ -76,11 +74,6 @@ namespace FluidSim.Core
link.Resolve(dtSub);
_timeOpenEnd += sw.Elapsed.TotalSeconds - t0;
- t0 = sw.Elapsed.TotalSeconds;
- foreach (var junc in _junctions)
- junc.Resolve(dtSub);
- _timeJunction += sw.Elapsed.TotalSeconds - t0;
-
t0 = sw.Elapsed.TotalSeconds;
foreach (var p in pipes)
p.SimulateSingleStep(dtSub);
diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs
index fbb1fae..9135ee3 100644
--- a/Core/SoundProcessor.cs
+++ b/Core/SoundProcessor.cs
@@ -3,53 +3,142 @@ using FluidSim.Core;
namespace FluidSim.Core
{
+ ///
+ /// Synthesises far‑field sound at a listener position from an open pipe end.
+ /// Three source mechanisms are combined:
+ /// 1. Monopole – time derivative of mass flow (dominant at low speed / high pulsation).
+ /// 2. Dipole – time derivative of momentum flux (shear‑layer / vortex shedding).
+ /// 3. Jet noise – Lighthill‑type turbulence mixing noise (scales with U^8).
+ ///
+ /// References:
+ /// • Lighthill, M.J. (1952) "On Sound Generated Aerodynamically".
+ /// • Dowling, A.P. & Williams, J.E.F. (1983) "Sound and Sources of Sound".
+ /// • Munjal, M.L. (2014) "Acoustics of Ducts and Mufflers", 2nd ed.
+ /// • Tam, C.K.W. & Auriault, L. (1999) "Jet Mixing Noise from Fine‑Scale Turbulence".
+ ///
public class SoundProcessor
{
private readonly double dt;
- private readonly double scaleFactor; // 1 / (4π r)
+ private readonly double c0; // ambient speed of sound (m/s)
+ private readonly double rho0; // ambient density (kg/m³)
+ private readonly double r; // listener distance (m)
+ private readonly double pipeArea; // cross‑sectional area of the pipe end (m²)
+
+ // ---------- monopole state ----------
+ private double flowLP;
+ private readonly double lpAlpha;
private double prevMassFlowOut;
private double smoothDMdt;
private readonly double alpha;
- // New: low‑pass the mass flow signal before derivative
- private double flowLP;
- private readonly double lpAlpha;
+ // ---------- dipole state ----------
+ private double prevMomentumFlux;
+ private double smoothDMomDt;
+ private readonly double dipAlpha;
+
+ // ---------- jet noise state ----------
+ private double jetNoiseSample; // previous random sample (for simple shaping)
+ private readonly double jetTau; // correlation time ≈ D / U_mean
public float Gain { get; set; } = 1.0f;
- public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0)
+ ///
+ ///
+ /// Audio sample rate (Hz).
+ /// Distance from the pipe exit to the listener (m).
+ /// Internal diameter of the pipe (m).
+ public SoundProcessor(int sampleRate,
+ double listenerDistanceMeters = 1.0,
+ double pipeDiameterMeters = 0.0217) // ~3.7 cm² area
{
dt = 1.0 / sampleRate;
- scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters);
+ r = listenerDistanceMeters;
+ pipeArea = Math.PI * 0.25 * pipeDiameterMeters * pipeDiameterMeters;
- // Smoothing time constant for the derivative: 10 ms (much smoother)
- double tau = 0.005; // 10 ms
+ // Ambient air properties
+ c0 = 340.0;
+ rho0 = 1.225;
+
+ // ---- Monopole smoothing ----
+ double tau = 0.002; // 2 ms
alpha = Math.Exp(-dt / tau);
- // Low‑pass time constant for the mass flow: 5 ms (kneecap high‑freq directly)
- double tauLP = 0.005;
+ double tauLP = 0.005; // 5 ms low‑pass on mass flow
lpAlpha = Math.Exp(-dt / tauLP);
+
+ // ---- Dipole smoothing ----
+ double tauDip = 0.003; // 3 ms
+ dipAlpha = Math.Exp(-dt / tauDip);
+
+ // ---- Jet noise correlation time ----
+ jetTau = Math.Max(0.0005, pipeDiameterMeters / 50.0); // D / U_ref, floor at 0.5 ms
}
+ ///
+ /// Process one sample. The OpenEndLink provides the instantaneous
+ /// exit‑plane mass flow, density, velocity, and pressure.
+ ///
public float Process(OpenEndLink openEnd)
{
- double flowOut = openEnd.LastMassFlowRate;
+ double flowOut = openEnd.LastMassFlowRate; // kg/s, positive = leaving pipe
+ double rhoExit = openEnd.LastFaceDensity; // kg/m³ at exit
+ double uExit = openEnd.LastFaceVelocity; // m/s (axial, positive = leaving)
+ double pExit = openEnd.LastFacePressure; // Pa
- // Low‑pass the mass flow signal
+ // ============================================================
+ // 1. MONOPOLE – due to unsteady mass addition (Lighthill 1952)
+ // Far‑field pressure: p'(r,t) = (1 / 4πr c0) · dṁ/dt
+ // ============================================================
flowLP = lpAlpha * flowLP + (1.0 - lpAlpha) * flowOut;
-
- // Derivative of the smoothed mass flow
double rawDerivative = (flowLP - prevMassFlowOut) / dt;
prevMassFlowOut = flowLP;
-
- // Smooth the derivative
smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative;
+ double pMono = smoothDMdt / (4.0 * Math.PI * r * c0);
- // Far‑field monopole pressure
- double pressure = smoothDMdt * scaleFactor * Gain;
+ // ============================================================
+ // 2. DIPOLE – due to unsteady momentum flux at the exit plane
+ // Momentum flux: F(t) = ṁ(t) · u(t) = ρ·A·u²
+ // Far‑field (compact, low M): p'(r,θ,t) ≈ (cosθ / 4πr c0) · dF/dt
+ // For on‑axis listener (θ = 0): p'(r,t) ≈ (1 / 4πr c0) · dF/dt
+ // We also include a U⁶ scaling factor relative to a reference velocity.
+ // ============================================================
+ double momentumFlux = Math.Abs(flowOut) * Math.Abs(uExit); // N
+ double rawMomDeriv = (momentumFlux - prevMomentumFlux) / dt;
+ prevMomentumFlux = momentumFlux;
+ smoothDMomDt = dipAlpha * smoothDMomDt + (1.0 - dipAlpha) * rawMomDeriv;
+ double pDipole = smoothDMomDt / (4.0 * Math.PI * r * c0);
- // Soft clip to ±1 (should rarely trigger now)
- return (float)pressure;
+ // Dipole efficiency factor: ∝ (U / c0)³ (since Idipole ∝ U⁶, pdipole ∝ U³)
+ double Mach = Math.Abs(uExit) / c0;
+ double dipoleEfficiency = Math.Pow(Mach, 3.0);
+ pDipole *= dipoleEfficiency;
+
+ // ============================================================
+ // 3. JET NOISE – Lighthill U⁸ mixing noise, band‑pass shaped
+ // rms pressure: p'_jet ~ ρ0 · A / r · U⁴ / c0²
+ // Model as broadband noise with amplitude ∝ U⁴.
+ // A simple first‑order low‑pass filter shapes the spectrum
+ // (cut‑off ≈ Strouhal frequency f ≈ 0.2 · U / D).
+ // ============================================================
+ double Uref = Math.Max(1.0, Math.Abs(uExit)); // avoid division by zero
+ double jetAmplitude = rho0 * pipeArea / r * Math.Pow(Uref / c0, 4.0);
+
+ // Correlation time (sample‑and‑hold style random walk)
+ double alphaJet = Math.Exp(-dt / jetTau);
+ // Generate a new random target each step, filter with alphaJet
+ double randomTarget = (new Random().NextDouble() * 2.0 - 1.0);
+ jetNoiseSample = alphaJet * jetNoiseSample + (1.0 - alphaJet) * randomTarget;
+ double pJet = jetAmplitude * jetNoiseSample;
+
+ // ============================================================
+ // Combine contributions (monopole is primary; dipole & jet are
+ // weighted down for realistic mix). Weights can be tuned per engine.
+ // ============================================================
+ double pressure = (3000.0 * pMono) + (0.01 * pDipole) + (0 * pJet);
+ pressure *= Gain;
+
+ // Soft‑clip to ±1
+ return (float)Math.Tanh(pressure);
}
}
}
\ No newline at end of file
diff --git a/Interfaces/Port.cs b/Interfaces/Port.cs
index 0b5a488..413c6b4 100644
--- a/Interfaces/Port.cs
+++ b/Interfaces/Port.cs
@@ -1,21 +1,14 @@
namespace FluidSim.Interfaces
{
- ///
- /// A fluid port that belongs to a component. The solver writes mass flow,
- /// enthalpy, etc. here, and reads pressure, density, etc.
- ///
public class Port
{
- // Set by the solver during coupling resolution
public double MassFlowRate; // kg/s, positive INTO the component that owns this port
- public double SpecificEnthalpy; // J/kg, enthalpy of the fluid crossing the port (inflow direction)
-
- // Set by the owning component after integration to reflect its current state
+ public double SpecificEnthalpy; // J/kg
public double Pressure; // Pa
public double Density; // kg/m³
public double Temperature; // K
+ public double AirFraction; // mass fraction of air (0 = exhaust, 1 = air)
- // Link back to owner (optional, but handy for debugging)
public object? Owner { get; set; }
public Port()
@@ -25,6 +18,7 @@
Pressure = 101325.0;
Density = 1.225;
Temperature = 300.0;
+ AirFraction = 1.0; // default fresh air
}
}
}
\ No newline at end of file
diff --git a/Program.cs b/Program.cs
index 11cd03d..c18f6ee 100644
--- a/Program.cs
+++ b/Program.cs
@@ -40,9 +40,7 @@ public class Program
// Throttle control
private static double _throttleTarget = 1.0; // 0‑1, set by arrow keys
private static double _throttleCurrent = 0.0; // actual current fraction (lerped)
- private const double ThrottleLerpRate = 5.0; // times per second (speed of movement)
- private const double ThrottleMinArea = 0.0000000000001; // 3.7e-5 m² ≈ 0.37 cm² (10% of pipe)
- private const double ThrottleMaxArea = 0.00000000001; // 3.7 cm² (full open)
+ private const double ThrottleLerpRate = 10.0; // times per second (speed of movement)
private static bool _wKeyHeld = false;
private static double _lastThrottleUpdateTime;
@@ -94,8 +92,7 @@ public class Program
_throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
}
- double actualArea = ThrottleMinArea + (ThrottleMaxArea - ThrottleMinArea) * _throttleCurrent;
- _scenario.ThrottleArea = actualArea;
+ _scenario.Throttle = _throttleCurrent;
// ---- Drawing ----
if (now - lastDrawTime >= 1.0 / DrawFrequency)
diff --git a/Scenarios/Scenario.cs b/Scenarios/Scenario.cs
index 0177104..f38e239 100644
--- a/Scenarios/Scenario.cs
+++ b/Scenarios/Scenario.cs
@@ -14,7 +14,30 @@ namespace FluidSim.Tests
protected const double AmbientPressure = 101325.0;
protected const double AmbientTemperature = 300.0;
- // ---------- Color helper ----------
+ // ---------- Color from pressure (volumes) ----------
+ protected Color PressureColor(double pressurePa)
+ {
+ double bar = pressurePa / 1e5; // convert to bar for easier mapping
+ byte r, g, b;
+
+ if (bar < 1.0) // vacuum → blue to green
+ {
+ double factor = Math.Clamp(bar, 0.0, 1.0);
+ r = 0;
+ g = (byte)(255 * factor);
+ b = (byte)(255 * (1.0 - factor));
+ }
+ else // above ambient → green to red
+ {
+ double factor = Math.Min((bar - 1.0) / 9.0, 1.0); // 1→10 bar maps to 0→1
+ r = (byte)(255 * factor);
+ g = (byte)(255 * (1.0 - factor));
+ b = 0;
+ }
+ return new Color(r, g, b);
+ }
+
+ // ---------- Color from temperature (pipes) ----------
protected Color TemperatureColor(double temperature)
{
double t = Math.Clamp(temperature, 0.0, 2000.0);
@@ -42,7 +65,7 @@ namespace FluidSim.Tests
{
var rect = new RectangleShape(new Vector2f(width, height))
{
- FillColor = TemperatureColor(volume.Temperature),
+ FillColor = PressureColor(volume.Pressure), // ← pressure‑based
Position = new Vector2f(centerX - width / 2f, topY)
};
target.Draw(rect);
@@ -60,7 +83,7 @@ namespace FluidSim.Tests
protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
float centerX, float topY, float width, float maxHeight)
{
- double fraction = cylinder.PistonFraction; // 0 = TDC, 1 = BDC
+ double fraction = cylinder.PistonFraction;
float currentHeight = (float)(maxHeight * fraction);
// Walls
@@ -69,10 +92,10 @@ namespace FluidSim.Tests
wall.Position = new Vector2f(centerX - width / 2f, topY);
target.Draw(wall);
- // Gas
+ // Gas – colored by pressure now
float gasTop = topY;
var gasRect = new RectangleShape(new Vector2f(width, currentHeight));
- gasRect.FillColor = TemperatureColor(cylinder.Temperature);
+ gasRect.FillColor = PressureColor(cylinder.Pressure); // ← pressure‑based
gasRect.Position = new Vector2f(centerX - width / 2f, gasTop);
target.Draw(gasRect);
@@ -95,7 +118,7 @@ namespace FluidSim.Tests
target.Draw(exhaustValve);
}
- // ---------- Draw a pipe ----------
+ // ---------- Draw a pipe (unchanged) ----------
protected void DrawPipe(RenderWindow target, Pipe1D pipe, float pipeCenterY, float pipeStartX, float pipeEndX)
{
int n = pipe.CellCount;
@@ -141,7 +164,7 @@ namespace FluidSim.Tests
{
float x = centers[i];
float r = radii[i];
- Color col = TemperatureColor(temperatures[i]);
+ Color col = TemperatureColor(temperatures[i]); // pipes still use temperature
stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col);
stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col);
diff --git a/Scenarios/TestScenario.cs b/Scenarios/TestScenario.cs
index 4454a3e..e7527c5 100644
--- a/Scenarios/TestScenario.cs
+++ b/Scenarios/TestScenario.cs
@@ -3,6 +3,7 @@ using SFML.Graphics;
using SFML.System;
using FluidSim.Components;
using FluidSim.Core;
+using FluidSim.Utils;
namespace FluidSim.Tests
{
@@ -10,67 +11,81 @@ namespace FluidSim.Tests
{
// Engine
private Cylinder cylinder;
+ private Crankshaft crankshaft;
// Intake side
- private Pipe1D intakePipeBeforeThrottle; // pipe from ambient to plenum
- private Volume0D intakePlenum; // plenum (100 mL)
- private Pipe1D intakeRunner; // pipe from plenum to cylinder
+ private Pipe1D intakePipeBeforeThrottle;
+ private Volume0D intakePlenum; // 5 mL
+ private Pipe1D intakeRunner;
// Exhaust side
private Pipe1D exhaustPipe;
// Links
- private OpenEndLink intakeOpenEnd; // ambient → left end of first pipe
- private OrificeLink throttleOrifice; // first pipe right end → plenum inlet (variable area)
- private OrificeLink plenumToRunner; // plenum outlet → runner left end (fixed area)
- private OrificeLink intakeValve; // runner right end → cylinder intake port
-
+ private OpenEndLink intakeOpenEnd;
+ private OrificeLink throttleOrifice;
+ private OrificeLink plenumToRunner;
+ private OrificeLink intakeValve;
private OrificeLink exhaustValve;
private OpenEndLink exhaustOpenEnd;
private Solver solver;
- private SoundProcessor soundProcessor;
+ private SoundProcessor exhaustSoundProcessor;
+ private SoundProcessor intakeSoundProcessor;
+ private OutdoorExhaustReverb reverb;
private double dt;
private int stepCount;
- public double ThrottleArea { get; set; } = 0.0; // controlled externally
+ // ---------- Throttle control ----------
+ public double Throttle { get; set; } = 0.0;
+ public double MaxThrottleArea { get; set; } = 6 * Units.cm2; // 2 cm²
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
- soundProcessor = new SoundProcessor(sampleRate, 1) { Gain = 1f };
solver = new Solver();
solver.SetTimeStep(dt);
solver.CflTarget = 0.9;
- // ---- Cylinder (no valve overlap to avoid backflow) ----
- double bore = 0.056, stroke = 0.050, conRod = 0.110, compRatio = 10.0;
+ // ---- Crankshaft (external, passed to cylinder) ----
+ crankshaft = new Crankshaft(1000);
+ crankshaft.Inertia = 0.05;
+ crankshaft.FrictionConstant = 2;
+ crankshaft.FrictionViscous = 0.05;
+
+ // ---- Cylinder ----
+ double bore = 0.056, stroke = 0.057, conRod = 0.110, compRatio = 9.2;
double ivo = 370.0, ivc = 580.0, evo = 120.0, evc = 350.0;
- cylinder = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, 1000)
+ cylinder = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft)
{
- MaxIntakeArea = 0.00037,
- MaxExhaustArea = 0.00037,
+ MaxIntakeArea = 3.7 * Units.cm2,
+ MaxExhaustArea = 3.7 * Units.cm2,
};
solver.AddComponent(cylinder);
- double pipeArea = 0.00037; // 3.7 cm²
+ double pipeDiameter = 2 * Units.cm;
+ double pipeArea = Units.AreaFromDiameter(pipeDiameter);
+
+ exhaustSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.05f };
+ intakeSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.05f };
+ reverb = new OutdoorExhaustReverb(sampleRate);
// ---- Pipes ----
- intakePipeBeforeThrottle = new Pipe1D(0.15, pipeArea, 5); // short pipe before throttle
- intakeRunner = new Pipe1D(0.10, pipeArea, 5); // runner after plenum
+ intakePipeBeforeThrottle = new Pipe1D(0.15, pipeArea, 5);
+ intakeRunner = new Pipe1D(0.1, pipeArea, 5);
exhaustPipe = new Pipe1D(1.00, pipeArea, 80);
solver.AddComponent(intakePipeBeforeThrottle);
solver.AddComponent(intakeRunner);
solver.AddComponent(exhaustPipe);
- // ---- Plenum (100 mL) ----
- intakePlenum = new Volume0D(0.0001, 101325.0, 300.0); // 0.0001 m³
- var plenumInlet = intakePlenum.CreatePort(); // from throttle
- var plenumOutlet = intakePlenum.CreatePort(); // to runner
+ // ---- Plenum (5 mL) ----
+ intakePlenum = new Volume0D(5 * Units.mL, 101325.0, 300.0);
+ var plenumInlet = intakePlenum.CreatePort();
+ var plenumOutlet = intakePlenum.CreatePort();
solver.AddComponent(intakePlenum);
- // ---- Intake open end (ambient → left end of first pipe) ----
+ // ---- Intake open end ----
intakeOpenEnd = new OpenEndLink(intakePipeBeforeThrottle, isLeftEnd: true)
{
AmbientPressure = 101325.0,
@@ -78,16 +93,16 @@ namespace FluidSim.Tests
};
solver.AddOpenEndLink(intakeOpenEnd);
- // ---- Throttle orifice (first pipe right end → plenum inlet) ----
+ // ---- Throttle orifice (variable area) ----
throttleOrifice = new OrificeLink(plenumInlet, intakePipeBeforeThrottle, isPipeLeftEnd: false,
- areaProvider: () => ThrottleArea)
+ areaProvider: () => MaxThrottleArea * Math.Clamp(Throttle, 0.001, 1))
{
- DischargeCoefficient = 0.1, // realistic throttle Cd
+ DischargeCoefficient = 0.1,
UseInertance = false
};
solver.AddOrificeLink(throttleOrifice);
- // ---- Plenum → runner (fixed area = pipe area) ----
+ // ---- Plenum to runner (fixed area) ----
plenumToRunner = new OrificeLink(plenumOutlet, intakeRunner, isPipeLeftEnd: true,
areaProvider: () => pipeArea)
{
@@ -96,7 +111,7 @@ namespace FluidSim.Tests
};
solver.AddOrificeLink(plenumToRunner);
- // ---- Intake valve (runner right end → cylinder intake port) ----
+ // ---- Intake valve ----
intakeValve = new OrificeLink(cylinder.IntakePort, intakeRunner, isPipeLeftEnd: false,
areaProvider: () => cylinder.IntakeValveArea)
{
@@ -130,16 +145,12 @@ namespace FluidSim.Tests
public override float Process()
{
- // 1. Advance crankshaft & pre‑step
cylinder.Crankshaft.Step(dt);
cylinder.PreStep(dt);
-
- // 2. Run solver
solver.Step();
stepCount++;
- // 3. Log every 200 steps
- if (stepCount % 200 == 0)
+ if (stepCount % 10000 == 0)
{
double crankDeg = cylinder.Crankshaft.CrankAngle * 180.0 / Math.PI % 720.0;
double cylP = cylinder.Pressure / 1e5;
@@ -149,14 +160,17 @@ namespace FluidSim.Tests
double mdotE = exhaustValve.LastMassFlowRate;
double pipeR = exhaustPipe.GetCellPressure(exhaustPipe.CellCount - 1) / 1e5;
double plenumP = intakePlenum.Pressure / 1e5;
+ double actualArea = MaxThrottleArea * Throttle;
Console.WriteLine($"Step {stepCount}: Angle={crankDeg:F1}°, " +
$"CylP={cylP:F2} bar, T={cylT:F0} K, mass={cylMass:F1} mg, " +
$"mdotI={mdotI:E4} kg/s, mdotE={mdotE:E4} kg/s, PipeR={pipeR:F2} bar");
- Console.WriteLine($"Throttle area = {ThrottleArea * 1e6:F2} mm², Plenum P = {plenumP:F3} bar");
+ Console.WriteLine($"Throttle = {Throttle * 100:F0}% area = {actualArea * 1e6:F2} mm², Plenum P = {plenumP:F3} bar");
}
- return soundProcessor.Process(exhaustOpenEnd);
+ float exhaustDry = exhaustSoundProcessor.Process(exhaustOpenEnd);
+ float intakeDry = intakeSoundProcessor.Process(intakeOpenEnd);
+ return reverb.Process(intakeDry + exhaustDry);
}
public override void Draw(RenderWindow target)
@@ -164,22 +178,21 @@ namespace FluidSim.Tests
float winW = target.GetView().Size.X;
float winH = target.GetView().Size.Y;
- // Fixed vertical centres for intake and exhaust
float intakeY = winH / 2f - 40f;
float exhaustY = winH / 2f + 80f;
- // ---- 1. Open end (ambient air source) ----
+ // Open end marker
float openEndX = 40f;
var openEndMark = new CircleShape(5f) { FillColor = Color.Cyan };
openEndMark.Position = new Vector2f(openEndX - 5f, intakeY - 5f);
target.Draw(openEndMark);
- // ---- 2. First intake pipe (ambient → throttle) ----
+ // First intake pipe
float pipe1StartX = openEndX;
float pipe1EndX = pipe1StartX + 120f;
DrawPipe(target, intakePipeBeforeThrottle, intakeY, pipe1StartX, pipe1EndX);
- // ---- 3. Throttle (symbolic restriction) ----
+ // Throttle symbol
float throttleX = pipe1EndX + 5f;
var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
{
@@ -188,25 +201,25 @@ namespace FluidSim.Tests
};
target.Draw(throttleRect);
- // ---- 4. Plenum ----
+ // Plenum
float plenW = 60f, plenH = 80f;
float plenLeftX = throttleX + 10f;
float plenCenterX = plenLeftX + plenW / 2f;
float plenTopY = intakeY - plenH / 2f;
DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
- // ---- 5. Runner pipe (plenum → cylinder) ----
+ // Runner pipe
float runnerStartX = plenLeftX + plenW + 5f;
float runnerEndX = runnerStartX + 100f;
DrawPipe(target, intakeRunner, intakeY, runnerStartX, runnerEndX);
- // ---- 6. Cylinder ----
- float cylCX = runnerEndX + 50f; // center X
- float cylTopY = intakeY - 120f; // top of cylinder (so it sits above the pipe)
+ // Cylinder
+ float cylCX = runnerEndX + 50f;
+ float cylTopY = intakeY - 120f;
float cylW = 80f, cylMaxH = 240f;
DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
- // ---- 7. Exhaust pipe (cylinder → open end) ----
+ // Exhaust pipe
float exhStartX = cylCX + cylW / 2f + 20f;
float exhEndX = winW - 60f;
DrawPipe(target, exhaustPipe, exhaustY, exhStartX, exhEndX);