From f79cf6b7ebba2359589241567ed240e5b4cf4b09 Mon Sep 17 00:00:00 2001 From: grillkol Date: Thu, 7 May 2026 13:28:41 +0200 Subject: [PATCH] orifice confirmed working --- Components/Volume0D.cs | 33 +++++----- Core/OpenEndLink.cs | 83 +++++++++++++----------- Core/Solver.cs | 42 ++++-------- Program.cs | 2 +- Scenarios/TestScenario.cs | 133 ++++++++++++++++++++------------------ 5 files changed, 143 insertions(+), 150 deletions(-) diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs index b7f4b13..fde55c0 100644 --- a/Components/Volume0D.cs +++ b/Components/Volume0D.cs @@ -4,23 +4,17 @@ using FluidSim.Interfaces; namespace FluidSim.Components { - /// - /// Zero‑dimensional control volume with arbitrary number of ports. - /// Integrates mass and energy fluxes from all ports. - /// Safeguards keep a tiny amount of gas to avoid negative states. - /// public class Volume0D : IComponent { public List Ports { get; } = new List(); - public double Mass { get; private set; } - public double InternalEnergy { get; private set; } + public double Mass { get; set; } // made public setter + public double InternalEnergy { get; set; } // made public setter public double Volume { get; set; } public double Dvdt { get; set; } public double Gamma { get; set; } = 1.4; public double GasConstant { get; set; } = 287.0; - // Ambient pressure used for emergency refill – default 101325 Pa public double AmbientPressure { get; set; } = 101325.0; // Derived quantities @@ -42,11 +36,9 @@ namespace FluidSim.Components InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0); } - /// Add a new port and initialise it to the volume's current state. public Port CreatePort() { var port = new Port { Owner = this }; - // Set the port state immediately to avoid a mismatch before the first integration port.Pressure = Pressure; port.Density = Density; port.Temperature = Temperature; @@ -56,9 +48,18 @@ namespace FluidSim.Components } /// - /// Integrate over dt using the MassFlowRate and SpecificEnthalpy - /// that have been set on each port during the coupling resolution phase. + /// Set the pressure to a specific value while keeping the current temperature constant. + /// Updates Mass and InternalEnergy accordingly. /// + public void SetPressure(double pressure) + { + double V = Math.Max(Volume, 1e-12); + double currentT = Temperature; // current temperature before changes + double rho = pressure / (GasConstant * currentT); + Mass = rho * V; + InternalEnergy = pressure * V / (Gamma - 1.0); + } + public void UpdateState(double dt) { double totalMdot = 0.0; @@ -67,17 +68,15 @@ namespace FluidSim.Components foreach (var port in Ports) { totalMdot += port.MassFlowRate; - // mdot * h gives energy flow: positive mdot = inflow, negative = outflow totalEdot += port.MassFlowRate * port.SpecificEnthalpy; } double dm = totalMdot * dt; - double dE = totalEdot * dt - Pressure * Dvdt * dt; // piston work + double dE = totalEdot * dt - Pressure * Dvdt * dt; Mass += dm; InternalEnergy += dE; - // Safeguards: keep at least 1 µg of gas at a small pressure double V = Math.Max(Volume, 1e-12); if (Mass < 1e-9) { @@ -89,18 +88,16 @@ namespace FluidSim.Components InternalEnergy = AmbientPressure * V / (Gamma - 1.0); } - // Final non‑negative clamps (should not be needed after above) if (Mass < 0.0) Mass = 1e-9; if (InternalEnergy < 0.0) InternalEnergy = AmbientPressure * V / (Gamma - 1.0); - // Push updated state back to all ports double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy; foreach (var port in Ports) { port.Pressure = p; port.Density = rho; port.Temperature = T; - port.SpecificEnthalpy = h; // will be overwritten by couplings for inflow, but this is the default + port.SpecificEnthalpy = h; } } diff --git a/Core/OpenEndLink.cs b/Core/OpenEndLink.cs index a6efce3..ca40d4a 100644 --- a/Core/OpenEndLink.cs +++ b/Core/OpenEndLink.cs @@ -3,14 +3,19 @@ using FluidSim.Components; namespace FluidSim.Core { + /// + /// Characteristic open‑end boundary condition after Jones (1978). + /// For all subsonic flow (outflow and inflow), the ghost state is derived + /// 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. + /// public class OpenEndLink { public Pipe1D Pipe { get; } public bool IsLeftEnd { get; } public double AmbientPressure { get; set; } = 101325.0; public double Gamma { get; set; } = 1.4; - public double GasConstant { get; set; } = 287.0; - public double AmbientTemperature { get; set; } = 300.0; public double LastMassFlowRate { get; private set; } public double LastFaceDensity { get; private set; } @@ -33,61 +38,63 @@ namespace FluidSim.Core double gm1 = gamma - 1.0; double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12)); double pAmb = AmbientPressure; - double rhoAmb = pAmb / (GasConstant * AmbientTemperature); - double aAmb = Math.Sqrt(gamma * pAmb / rhoAmb); + + // Riemann invariants + double J_plus = uInt + 2.0 * cInt / gm1; + double J_minus = uInt - 2.0 * cInt / gm1; double rhoGhost, uGhost, pGhost; - // ----- Supersonic outflow: extrapolate interior ----- - bool supersonicOut = IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt); - if (supersonicOut) + // ---- 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); + + // Check for supersonic outflow: if the isentropic velocity exceeds the speed of sound, + // the flow is supersonic and we extrapolate the interior state. + bool supersonic = IsLeftEnd + ? (uInt <= -cInt) // left end: outflow is when u < -c + : (uInt >= cInt); // right end: outflow is when u > c + + // 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) + supersonic = uIso <= -cIso; + else + supersonic = uIso >= cIso; + } + + if (supersonic) + { + // Supersonic outflow – extrapolate interior rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; } else { - // Riemann invariants - double J_plus = uInt + 2.0 * cInt / gm1; - double J_minus = uInt - 2.0 * cInt / gm1; - - // Trial subsonic outflow ghost state - double s = pInt / Math.Pow(rhoInt, gamma); - double rhoOut = Math.Pow(pAmb / s, 1.0 / gamma); - double cOut = Math.Sqrt(gamma * pAmb / rhoOut); - double uOut = IsLeftEnd - ? (J_minus + 2.0 * cOut / gm1) - : (J_plus - 2.0 * cOut / gm1); - - bool outflowPossible = IsLeftEnd ? (uOut <= 0) : (uOut >= 0); - - if (outflowPossible) - { - // Subsonic outflow - pGhost = pAmb; - rhoGhost = rhoOut; - uGhost = uOut; - } - else - { - // Subsonic inflow (ambient reservoir model) - pGhost = pAmb; - rhoGhost = rhoAmb; - uGhost = IsLeftEnd - ? (J_minus + 2.0 * aAmb / gm1) - : (J_plus - 2.0 * aAmb / gm1); - } + // Subsonic flow – use the isentropic state + rhoGhost = rhoIso; + uGhost = uIso; + pGhost = pAmb; } + // Apply ghost to pipe if (IsLeftEnd) Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost); else Pipe.SetGhostRight(rhoGhost, uGhost, pGhost); + // Mass flow out of the pipe (positive = leaving) double area = Pipe.Area; double mdot = rhoGhost * uGhost * area; - if (IsLeftEnd) mdot = -mdot; + if (IsLeftEnd) mdot = -mdot; // left end: positive u is into pipe, so out is -u LastMassFlowRate = mdot; LastFaceDensity = rhoGhost; LastFaceVelocity = uGhost; diff --git a/Core/Solver.cs b/Core/Solver.cs index c3d348b..40ba7e9 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -6,10 +6,6 @@ using FluidSim.Interfaces; namespace FluidSim.Core { - /// - /// Top‑level solver that owns all components and couplings, - /// orchestrates sub‑stepping, and exposes states for audio. - /// public class Solver { private readonly List _components = new(); @@ -19,6 +15,9 @@ namespace FluidSim.Core private double _dt; + /// CFL target for sub‑stepping (0.3‑0.8). Lower values are safer for shocks. + public double CflTarget { get; set; } = 0.8; + public void SetTimeStep(double dt) => _dt = dt; public void AddComponent(IComponent component) => _components.Add(component); @@ -26,57 +25,38 @@ namespace FluidSim.Core public void AddJunction(Junction junction) => _junctions.Add(junction); public void AddOpenEndLink(OpenEndLink link) => _openEndLinks.Add(link); - // Convenience: first pipe’s port B mass flow (often the exhaust) - public double ExhaustMassFlow - { - get - { - var pipes = _components.OfType().ToList(); - if (pipes.Count > 0) - return Math.Abs(pipes[0].PortB.MassFlowRate); - return 0.0; - } - } - - /// - /// Advance the whole system by one global time step. - /// public void Step() { var pipes = _components.OfType().ToList(); if (pipes.Count == 0) return; - // 1. Determine sub‑step count (max CFL over all pipes) int nSub = 1; foreach (var p in pipes) - nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt)); + nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt, CflTarget)); double dtSub = _dt / nSub; - // 2. Sub‑step loop + const int maxSubSteps = 10000; + if (nSub > maxSubSteps) + { + Console.WriteLine($"Warning: required sub‑steps {nSub} exceeds limit. Simulation stopped."); + return; + } + for (int sub = 0; sub < nSub; sub++) { - // a) Resolve all orifice links (volume ↔ pipe) foreach (var link in _orificeLinks) link.Resolve(dtSub); - - // b) Resolve all open‑end links (pipe → atmosphere) foreach (var link in _openEndLinks) link.Resolve(dtSub); - - // c) Resolve all junctions (pipe ↔ pipe) foreach (var junc in _junctions) junc.Resolve(dtSub); - - // d) Advance all pipes foreach (var p in pipes) p.SimulateSingleStep(dtSub); } - // 3. Clear ghost flags foreach (var p in pipes) p.ClearGhostFlags(); - // 4. Integrate non‑pipe components (volumes, atmosphere, etc.) foreach (var comp in _components) comp.UpdateState(_dt); } diff --git a/Program.cs b/Program.cs index f4cd113..74e078e 100644 --- a/Program.cs +++ b/Program.cs @@ -14,7 +14,7 @@ public class Program private static Scenario scenario; // Speed control (existing + new throttle) - private static double desiredSpeed = 0.001; + private static double desiredSpeed = 0.01; //private static double desiredSpeed = 1.0; private static double currentSpeed = desiredSpeed; private const double MinSpeed = 0.0001; diff --git a/Scenarios/TestScenario.cs b/Scenarios/TestScenario.cs index 6178091..32cc958 100644 --- a/Scenarios/TestScenario.cs +++ b/Scenarios/TestScenario.cs @@ -3,60 +3,74 @@ using SFML.Graphics; using SFML.System; using FluidSim.Components; using FluidSim.Core; +using FluidSim.Utils; namespace FluidSim.Tests { public class TestScenario : Scenario { private Solver solver; + private Volume0D volume; private Pipe1D pipe; - private OrificeLink closedEnd; // left end – closed wall - private OpenEndLink openEndLink; // right end – atmosphere + private OrificeLink orifice; // volume → pipe left + private OpenEndLink openEnd; // pipe right → atmosphere private SoundProcessor soundProcessor; - private OutdoorExhaustReverb reverb; private int stepCount; - private double simTime; // elapsed simulation time (seconds) - private double pulseInterval = 0.1; // seconds between pulses - private double nextPulseTime; + + // Pressure reset scheduling + private double simTime; private double dt; + private double resetInterval = 0.2; // seconds between resets + private double nextResetTime; + private double targetPressure = 10 * Units.atm; + private double rampDuration = 0.002; // 2 ms ramp + private double rampStartTime; + private double rampStartPressure; // pressure when ramp started + private bool ramping; public override void Initialize(int sampleRate) { dt = 1.0 / sampleRate; soundProcessor = new SoundProcessor(sampleRate, 1); - soundProcessor.Gain = 10; - - reverb = new OutdoorExhaustReverb(sampleRate); + soundProcessor.Gain = 2.0f; // lower gain to avoid clipping solver = new Solver(); solver.SetTimeStep(dt); + solver.CflTarget = 0.4; - // Pipe: 2 m long, 1 cm² area, 200 cells - pipe = new Pipe1D(length: 2, area: 1e-4, cellCount: 100); + volume = new Volume0D(1e-3, targetPressure, 300.0); + solver.AddComponent(volume); + + pipe = new Pipe1D(1.0, 1e-4, 200); + pipe.EnergyRelaxationRate = 10; solver.AddComponent(pipe); - // Initially pipe at ambient conditions - pipe.SetUniformState(1.225, 0.0, 101325.0); + var volPort = volume.CreatePort(); + double orificeArea = 1e-5; + orifice = new OrificeLink(volPort, pipe, isPipeLeftEnd: true, + areaProvider: () => orificeArea) + { + DischargeCoefficient = 0.62, + UseInertance = true, + EffectiveLength = 0.001 + }; + solver.AddOrificeLink(orifice); - // Left end: closed wall (area = 0 → reflective) - closedEnd = new OrificeLink(null, pipe, isPipeLeftEnd: true, areaProvider: () => 0.0); - solver.AddOrificeLink(closedEnd); - - // Right end: open to atmosphere - openEndLink = new OpenEndLink(pipe, isLeftEnd: false) + openEnd = new OpenEndLink(pipe, isLeftEnd: false) { AmbientPressure = 101325.0, Gamma = 1.4 }; - solver.AddOpenEndLink(openEndLink); + solver.AddOpenEndLink(openEnd); stepCount = 0; simTime = 0.0; - nextPulseTime = pulseInterval; // first pulse at 100 ms + nextResetTime = resetInterval; + ramping = false; - Console.WriteLine("Pulse reflection test – closed left, open right"); - Console.WriteLine("Pulse injected every 100 ms at left end (cell 0)"); + Console.WriteLine("Pressure reset test with smooth ramp"); + Console.WriteLine($"Volume 1L, reset to {targetPressure} Pa every {resetInterval*1000} ms, ramp {rampDuration*1000} ms"); } public override float Process() @@ -65,59 +79,54 @@ namespace FluidSim.Tests stepCount++; simTime += dt; - // ---- Inject a pressure pulse at the closed end every 100 ms ---- - if (simTime >= nextPulseTime) + // ---- Smooth pressure ramp ---- + if (ramping) { - // Apply a Gaussian pulse to the first few cells - double ambientPressure = 101325.0; - double pulseAmplitude = 20 * ambientPressure; // 0.5 atm overpressure - double pulseWidth = 0.05; // m (spread over a few cells) - int n = pipe.CellCount; - double dx = 2.0 / n; - - // Only modify cells within 2*pulseWidth from the left end - int maxCell = Math.Min(5, n - 1); // at most the first 5 cells - for (int i = 0; i <= maxCell; i++) + double elapsed = simTime - rampStartTime; + if (elapsed >= rampDuration) { - double x = (i + 0.5) * dx; - double P = pulseAmplitude * Math.Exp(-x * x / (pulseWidth * pulseWidth)); - double currentP = pipe.GetCellPressure(i); - double newP = P; - // Update pressure, keeping density and velocity unchanged - // We recompute total energy accordingly - double rho = pipe.GetCellDensity(i); - double u = pipe.GetCellVelocity(i); - double e = newP / ((1.4 - 1.0) * rho); - double E = rho * e + 0.5 * rho * u * u; - pipe.SetCellState(i, rho, u, newP); + // Ramp finished, set exactly to target + volume.SetPressure(targetPressure); + ramping = false; + } + else + { + double frac = elapsed / rampDuration; + double currentTarget = rampStartPressure + (targetPressure - rampStartPressure) * frac; + volume.SetPressure(currentTarget); } - Console.WriteLine($"Pulse injected at t = {simTime:F3} s"); - nextPulseTime += pulseInterval; } - // Audio from open‑end mass flow - float sample = soundProcessor.Process(openEndLink); - - // Log every 200 steps - if (stepCount % 1000 == 0) + // ---- Trigger a new reset ---- + if (!ramping && simTime >= nextResetTime) { - int leftIdx = 0; - int midIdx = pipe.CellCount / 2; - int rightIdx = pipe.CellCount - 1; - double pL = pipe.GetCellPressure(leftIdx); - double pM = pipe.GetCellPressure(midIdx); - double pR = pipe.GetCellPressure(rightIdx); + rampStartPressure = volume.Pressure; // current pressure before reset + rampStartTime = simTime; + ramping = true; + nextResetTime += resetInterval; + } - Console.WriteLine($"Step {stepCount}: P_left={pL:F1} Pa, P_mid={pM:F1} Pa, P_right={pR:F1} Pa"); + // Log every 500 steps + if (stepCount % 500 == 0) + { + double volP = volume.Pressure; + double pipeL = pipe.GetCellPressure(0); + double pipeR = pipe.GetCellPressure(pipe.CellCount - 1); + double mdotOrif = orifice.LastMassFlowRate; + double mdotOpen = openEnd.LastMassFlowRate; + + Console.WriteLine($"Step {stepCount}: " + + $"VolP={volP:F1} Pa, PipeL={pipeL:F1}, PipeR={pipeR:F1}, " + + $"mdot_orif={mdotOrif:E4} kg/s, mdot_open={mdotOpen:E4} kg/s"); } if (double.IsNaN(pipe.GetCellPressure(0))) { - Console.WriteLine("NaN detected – stopping simulation."); + Console.WriteLine("NaN detected – stopping."); return 0f; } - return reverb.Process(sample); + return soundProcessor.Process(openEnd); } public override void Draw(RenderWindow target)