From b3230844b7d41d18a66dd3f692d1c41d2d5c505b Mon Sep 17 00:00:00 2001 From: max Date: Thu, 7 May 2026 21:48:37 +0200 Subject: [PATCH] Engine working --- Components/Atmosphere.cs | 2 +- Components/Cylinder.cs | 122 ++++++++++--------- Components/Pipe1D.cs | 129 ++++++++++++-------- Components/Volume0D.cs | 52 +++++--- Core/Junction.cs | 228 ----------------------------------- Core/OpenEndLink.cs | 52 +++++--- Core/OrificeLink.cs | 43 ++++--- Core/OutdoorExhaustReverb.cs | 2 +- Core/Solver.cs | 7 -- Core/SoundProcessor.cs | 129 +++++++++++++++++--- Interfaces/Port.cs | 12 +- Program.cs | 7 +- Scenarios/Scenario.cs | 37 ++++-- Scenarios/TestScenario.cs | 105 +++++++++------- 14 files changed, 441 insertions(+), 486 deletions(-) delete mode 100644 Core/Junction.cs 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);