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

Engine working

Browse Source
This commit is contained in:
max
2026-05-07 21:48:37 +02:00
parent 92d84eacfe
commit b3230844b7
14 changed files with 441 additions and 486 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Components/Atmosphere.cs
View File

@@ -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;
}
}
}

122
Components/Cylinder.cs
View File

@@ -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;
}
}
}

129
Components/Pipe1D.cs
View File

@@ -7,10 +7,10 @@ namespace FluidSim.Components
/// <summary>
/// 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.
/// </summary>
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)

52
Components/Volume0D.cs
View File

@@ -1,4 +1,4 @@
using System;
using System;
using System.Collections.Generic;
using FluidSim.Interfaces;
@@ -8,8 +8,9 @@ namespace FluidSim.Components
{
public List<Port> Ports { get; } = new List<Port>();
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;
}
/// <summary>
/// Set the pressure to a specific value while keeping the current temperature constant.
/// Updates Mass and InternalEnergy accordingly.
/// </summary>
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;
}
}

228
Core/Junction.cs
View File

@@ -1,228 +0,0 @@
using System;
using System.Collections.Generic;
using FluidSim.Components;
namespace FluidSim.Core
{
/// <summary>
/// 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.
/// </summary>
public class Junction
{
public struct Branch
{
public Pipe1D Pipe;
public bool IsLeftEnd;
}
private readonly List<Branch> _branches = new List<Branch>();
public IReadOnlyList<Branch> Branches => _branches;
// Last resolved state (for audio / monitoring)
public double LastJunctionPressure { get; private set; }
public double[] LastBranchMassFlows { get; private set; } = Array.Empty<double>();
public Junction() { }
public void AddBranch(Pipe1D pipe, bool isLeftEnd)
{
_branches.Add(new Branch { Pipe = pipe, IsLeftEnd = isLeftEnd });
}
/// <summary>
/// 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.
/// </summary>
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<double, double> 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);
}
}
/// <summary>Simple Brent's method root finder.</summary>
private static double BrentsMethod(Func<double, double> 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;
}
}
}

52
Core/OpenEndLink.cs
View File

@@ -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.
/// </summary>
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;

43
Core/OrificeLink.cs
View File

@@ -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;

2
Core/OutdoorExhaustReverb.cs
View File

@@ -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)

7
Core/Solver.cs
View File

@@ -11,7 +11,6 @@ namespace FluidSim.Core
{
private readonly List<IComponent> _components = new();
private readonly List<OrificeLink> _orificeLinks = new();
private readonly List<Junction> _junctions = new();
private readonly List<OpenEndLink> _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);

129
Core/SoundProcessor.cs
View File

@@ -3,53 +3,142 @@ using FluidSim.Core;
namespace FluidSim.Core
{
/// <summary>
/// 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".
/// </summary>
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)
/// <summary>
/// </summary>
/// <param name="sampleRate">Audio sample rate (Hz).</param>
/// <param name="listenerDistanceMeters">Distance from the pipe exit to the listener (m).</param>
/// <param name="pipeDiameterMeters">Internal diameter of the pipe (m).</param>
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
}
/// <summary>
/// Process one sample. The OpenEndLink provides the instantaneous
/// exit‑plane mass flow, density, velocity, and pressure.
/// </summary>
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);
}
}
}

12
Interfaces/Port.cs
View File

@@ -1,21 +1,14 @@
namespace FluidSim.Interfaces
{
/// <summary>
/// A fluid port that belongs to a component. The solver writes mass flow,
/// enthalpy, etc. here, and reads pressure, density, etc.
/// </summary>
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
}
}
}

7
Program.cs
View File

@@ -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)

37
Scenarios/Scenario.cs
View File

@@ -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);

105
Scenarios/TestScenario.cs
View File

@@ -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);