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)