diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
index a77b885..fd19cb3 100644
--- a/Components/Pipe1D.cs
+++ b/Components/Pipe1D.cs
@@ -3,24 +3,35 @@ using FluidSim.Interfaces;
namespace FluidSim.Components
{
+ public enum BoundaryType
+ {
+ VolumeCoupling,
+ OpenEnd,
+ ClosedEnd
+ }
+
public class Pipe1D
{
public Port PortA { get; }
public Port PortB { get; }
public double Area => _area;
- private int _n; // number of cells
+ private int _n;
private double _dx, _dt, _gamma, _area;
private double[] _rho, _rhou, _E;
- // Volume boundary states, constant during sub‑steps
private double _rhoLeft, _pLeft;
private double _rhoRight, _pRight;
private bool _leftBCSet, _rightBCSet;
- // CFL control
+ private BoundaryType _leftBCType = BoundaryType.VolumeCoupling;
+ private BoundaryType _rightBCType = BoundaryType.VolumeCoupling;
+
+ private double _leftAmbientPressure = 101325.0;
+ private double _rightAmbientPressure = 101325.0;
+
private const double CflTarget = 0.8;
- private const double ReferenceSoundSpeed = 340.0; // m/s, standard air
+ private const double ReferenceSoundSpeed = 340.0;
public double FrictionFactor { get; set; } = 0.02;
@@ -29,28 +40,29 @@ namespace FluidSim.Components
public double GetCellPressure(int i) => Pressure(i);
public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
- ///
- /// Creates a 1D pipe.
- /// Cell count is automatically determined to satisfy CFL in still air.
- ///
- /// Pipe length in metres.
- /// Cross‑sectional area in m².
- /// Global simulation sample rate (Hz).
- public Pipe1D(double length, double area, int sampleRate)
- {
- // Desired spatial step to keep CFL ≤ target for still air
- double dtGlobal = 1.0 / sampleRate;
- double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
+ public BoundaryType LeftBCType => _leftBCType;
+ public BoundaryType RightBCType => _rightBCType;
- // Number of cells must be at least 2; try to hit dxTarget
- int nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
- // Ensure we don't accidentally overshoot dxTarget by more than a factor
- while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
- nCells++;
+ public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
+ {
+ double dtGlobal = 1.0 / sampleRate;
+ int nCells;
+
+ if (forcedCellCount > 1)
+ {
+ nCells = forcedCellCount;
+ }
+ else
+ {
+ double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
+ nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
+ while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
+ nCells++;
+ }
_n = nCells;
_dx = length / _n;
- _dt = dtGlobal; // global (audio) time step
+ _dt = dtGlobal;
_area = area;
_gamma = 1.4;
@@ -62,6 +74,11 @@ namespace FluidSim.Components
PortB = new Port();
}
+ public void SetLeftBoundaryType(BoundaryType type) => _leftBCType = type;
+ public void SetRightBoundaryType(BoundaryType type) => _rightBCType = type;
+ public void SetLeftAmbientPressure(double p) => _leftAmbientPressure = p;
+ public void SetRightAmbientPressure(double p) => _rightAmbientPressure = p;
+
public void SetUniformState(double rho, double u, double p)
{
double e = p / ((_gamma - 1) * rho);
@@ -74,10 +91,14 @@ namespace FluidSim.Components
}
}
- public double GetLeftPressure() => Pressure(0);
- public double GetRightPressure() => Pressure(_n - 1);
- public double GetLeftDensity() => _rho[0];
- public double GetRightDensity() => _rho[_n - 1];
+ public void SetCellState(int i, double rho, double u, double p)
+ {
+ if (i < 0 || i >= _n) return;
+ _rho[i] = rho;
+ _rhou[i] = rho * u;
+ double e = p / ((_gamma - 1) * rho);
+ _E[i] = rho * e + 0.5 * rho * u * u;
+ }
public void SetLeftVolumeState(double rhoVol, double pVol)
{
@@ -93,6 +114,165 @@ namespace FluidSim.Components
_rightBCSet = true;
}
+ public void ClearBC() => _leftBCSet = _rightBCSet = false;
+
+ public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
+ {
+ double maxW = 0.0;
+ for (int i = 0; i < _n; i++)
+ {
+ double rho = _rho[i];
+ double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
+ double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
+ double local = u + c;
+ if (local > maxW) maxW = local;
+ }
+ maxW = Math.Max(maxW, 1e-8);
+ return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
+ }
+
+ public void SimulateSingleStep(double dtSub)
+ {
+ int n = _n;
+ double[] Fm = new double[n + 1];
+ double[] Fp = new double[n + 1];
+ double[] Fe = new double[n + 1];
+
+ // Left boundary (face 0)
+ switch (_leftBCType)
+ {
+ case BoundaryType.VolumeCoupling:
+ if (_leftBCSet)
+ {
+ HLLCFlux(_rhoLeft, 0.0, _pLeft,
+ _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), Pressure(0),
+ out Fm[0], out Fp[0], out Fe[0]);
+ }
+ else
+ {
+ Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0;
+ }
+ break;
+
+ case BoundaryType.OpenEnd:
+ {
+ double rhoR = _rho[0];
+ double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
+ double pR = Pressure(0);
+ HLLCFlux(rhoR, uR, _leftAmbientPressure,
+ rhoR, uR, pR,
+ out Fm[0], out Fp[0], out Fe[0]);
+ }
+ break;
+
+ case BoundaryType.ClosedEnd:
+ {
+ double rhoR = _rho[0];
+ double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
+ double pR = Pressure(0);
+ HLLCFlux(rhoR, -uR, pR,
+ rhoR, uR, pR,
+ out Fm[0], out Fp[0], out Fe[0]);
+ }
+ break;
+ }
+
+ // Internal faces
+ for (int i = 0; i < n - 1; i++)
+ {
+ double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
+ double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
+ HLLCFlux(_rho[i], uL, Pressure(i),
+ _rho[i + 1], uR, Pressure(i + 1),
+ out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
+ }
+
+ // Right boundary (face n)
+ switch (_rightBCType)
+ {
+ case BoundaryType.VolumeCoupling:
+ if (_rightBCSet)
+ {
+ double rhoL = _rho[n - 1];
+ double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
+ double pL = Pressure(n - 1);
+ HLLCFlux(rhoL, uL, pL,
+ _rhoRight, 0.0, _pRight,
+ out Fm[n], out Fp[n], out Fe[n]);
+ }
+ else
+ {
+ Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0;
+ }
+ break;
+
+ case BoundaryType.OpenEnd:
+ {
+ double rhoL = _rho[n - 1];
+ double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
+ double pL = Pressure(n - 1);
+ HLLCFlux(rhoL, uL, pL,
+ rhoL, uL, _rightAmbientPressure,
+ out Fm[n], out Fp[n], out Fe[n]);
+ }
+ break;
+
+ case BoundaryType.ClosedEnd:
+ {
+ double rhoL = _rho[n - 1];
+ double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
+ double pL = Pressure(n - 1);
+ HLLCFlux(rhoL, uL, pL,
+ rhoL, -uL, pL,
+ out Fm[n], out Fp[n], out Fe[n]);
+ }
+ break;
+ }
+
+ // Cell update
+ for (int i = 0; i < n; i++)
+ {
+ double dM = (Fm[i + 1] - Fm[i]) / _dx;
+ double dP = (Fp[i + 1] - Fp[i]) / _dx;
+ double dE = (Fe[i + 1] - Fe[i]) / _dx;
+
+ _rho[i] -= dtSub * dM;
+ _rhou[i] -= dtSub * dP;
+ _E[i] -= dtSub * dE;
+
+ if (_rho[i] < 1e-12) _rho[i] = 1e-12;
+
+ double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
+ double pMin = 100.0;
+ double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
+ if (_E[i] < eMin) _E[i] = eMin;
+ }
+
+ // Port quantities (only meaningful for volume coupled ends)
+ double mdotA_sub = _leftBCType == BoundaryType.VolumeCoupling && _leftBCSet ? Fm[0] * _area : 0.0;
+ double mdotB_sub = _rightBCType == BoundaryType.VolumeCoupling && _rightBCSet ? -Fm[n] * _area : 0.0;
+
+ PortA.MassFlowRate = mdotA_sub;
+ PortB.MassFlowRate = mdotB_sub;
+ PortA.Pressure = Pressure(0);
+ PortB.Pressure = Pressure(_n - 1);
+ PortA.Density = _rho[0];
+ PortB.Density = _rho[_n - 1];
+
+ if (_leftBCType == BoundaryType.VolumeCoupling && _leftBCSet)
+ {
+ PortA.SpecificEnthalpy = mdotA_sub < 0
+ ? GetCellTotalSpecificEnthalpy(0)
+ : 0.0;
+ }
+ if (_rightBCType == BoundaryType.VolumeCoupling && _rightBCSet)
+ {
+ PortB.SpecificEnthalpy = mdotB_sub < 0
+ ? GetCellTotalSpecificEnthalpy(_n - 1)
+ : 0.0;
+ }
+ }
+
private double GetCellTotalSpecificEnthalpy(int i)
{
double rho = Math.Max(_rho[i], 1e-12);
@@ -102,146 +282,6 @@ namespace FluidSim.Components
return h + 0.5 * u * u;
}
- ///
- /// Advance the pipe over one global time step using sub‑stepping.
- /// Must be called once per global simulation cycle.
- ///
- public void Simulate()
- {
- int n = _n;
-
- // --- Determine maximum wave speed in the pipe ---
- double maxWaveSpeed = 0.0;
- for (int i = 0; i < n; i++)
- {
- double rho = Math.Max(_rho[i], 1e-12);
- double u = Math.Abs(_rhou[i] / rho);
- double c = Math.Sqrt(_gamma * Pressure(i) / rho);
- double local = u + c;
- if (local > maxWaveSpeed) maxWaveSpeed = local;
- }
-
- if (maxWaveSpeed < 1e-8) maxWaveSpeed = 1e-8;
-
- int nSub = Math.Max(1, (int)Math.Ceiling(_dt * maxWaveSpeed / (CflTarget * _dx)));
- double dtSub = _dt / nSub;
-
- // Accumulators for net mass flows
- double sumMdotA = 0.0, sumMdotB = 0.0;
-
- // Accumulators for fluid that ENTERS the volumes (pipe → volume)
- double massInA = 0.0, energyInA = 0.0;
- double massInB = 0.0, energyInB = 0.0;
-
- for (int step = 0; step < nSub; step++)
- {
- double[] Fm = new double[n + 1];
- double[] Fp = new double[n + 1];
- double[] Fe = new double[n + 1];
-
- // Left boundary (face 0)
- if (_leftBCSet)
- {
- HLLCFlux(_rhoLeft, 0.0, _pLeft,
- _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), Pressure(0),
- out Fm[0], out Fp[0], out Fe[0]);
- }
- else
- {
- Fm[0] = 0;
- Fp[0] = Pressure(0);
- Fe[0] = 0;
- }
-
- // Internal faces
- for (int i = 0; i < n - 1; i++)
- {
- double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
- double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
- HLLCFlux(_rho[i], uL, Pressure(i),
- _rho[i + 1], uR, Pressure(i + 1),
- out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
- }
-
- // Right boundary (face n)
- if (_rightBCSet)
- {
- double rhoL = _rho[n - 1];
- double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
- double pL = Pressure(n - 1);
- HLLCFlux(rhoL, uL, pL,
- _rhoRight, 0.0, _pRight,
- out Fm[n], out Fp[n], out Fe[n]);
- }
- else
- {
- Fm[n] = 0;
- Fp[n] = Pressure(n - 1);
- Fe[n] = 0;
- }
-
- // Cell update
- for (int i = 0; i < n; i++)
- {
- double dM = (Fm[i + 1] - Fm[i]) / _dx;
- double dP = (Fp[i + 1] - Fp[i]) / _dx;
- double dE = (Fe[i + 1] - Fe[i]) / _dx;
-
- _rho[i] -= dtSub * dM;
- _rhou[i] -= dtSub * dP;
- _E[i] -= dtSub * dE;
-
- if (_rho[i] < 1e-12) _rho[i] = 1e-12;
- double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
- if (_E[i] < kinetic) _E[i] = kinetic;
- }
-
- // Sub‑step mass flow rates (kg/s)
- double mdotA_sub = _leftBCSet ? Fm[0] * _area : 0.0; // >0 = into pipe
- double mdotB_sub = _rightBCSet ? -Fm[n] * _area : 0.0; // >0 = into pipe from right
-
- sumMdotA += mdotA_sub;
- sumMdotB += mdotB_sub;
-
- // Flow FROM pipe INTO volume A: mdotA_sub < 0
- if (mdotA_sub < 0 && _leftBCSet)
- {
- double massRate = -mdotA_sub; // kg/s entering volume A
- double h = GetCellTotalSpecificEnthalpy(0);
- massInA += massRate * dtSub;
- energyInA += massRate * dtSub * h;
- }
-
- // Flow FROM pipe INTO volume B: mdotB_sub < 0 (because
- // mdotB_sub = -Fm[n], and Fm[n] > 0 is flow to the right)
- if (mdotB_sub < 0 && _rightBCSet)
- {
- double massRate = -mdotB_sub; // kg/s entering volume B
- double h = GetCellTotalSpecificEnthalpy(_n - 1);
- massInB += massRate * dtSub;
- energyInB += massRate * dtSub * h;
- }
- }
-
- // Averaged net mass flows (sign: positive = into pipe)
- PortA.MassFlowRate = sumMdotA / nSub;
- PortB.MassFlowRate = sumMdotB / nSub;
-
- // Assign enthalpy ONLY for the fluid that physically entered the volume
- if (massInA > 1e-12)
- PortA.SpecificEnthalpy = energyInA / massInA;
-
- if (massInB > 1e-12)
- PortB.SpecificEnthalpy = energyInB / massInB;
-
- // If no inflow occurred, leave the port’s enthalpy unchanged.
- // (It will be set to the volume’s static enthalpy by PushStateToPort
- // or overwritten by TransferPipeToVolume if flow reverses later.)
-
- _leftBCSet = _rightBCSet = false;
- }
-
- // Pressure and HLLC flux unchanged
private double Pressure(int i) =>
(_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
@@ -279,5 +319,12 @@ namespace FluidSim.Components
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
}
}
+
+ public double GetPressureAtFraction(double fraction)
+ {
+ int i = (int)(fraction * (_n - 1));
+ i = Math.Clamp(i, 0, _n - 1);
+ return Pressure(i);
+ }
}
}
\ No newline at end of file
diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs
index b5949a2..2ff3940 100644
--- a/Components/Volume0D.cs
+++ b/Components/Volume0D.cs
@@ -50,13 +50,26 @@ namespace FluidSim.Components
Port.SpecificEnthalpy = SpecificEnthalpy;
}
+ // Original integrate (uses the constructor’s sample rate)
public void Integrate()
+ {
+ Integrate(_dt);
+ }
+
+ public void SetPressure(double newPressure)
+ {
+ InternalEnergy = newPressure * Volume / (Gamma - 1.0);
+ // Mass stays the same, so density is unchanged
+ }
+
+ // New overload: integrate with a custom time step (for sub‑steps)
+ public void Integrate(double dtOverride)
{
double mdot = Port.MassFlowRate;
double h_in = Port.SpecificEnthalpy;
- double dm = mdot * _dt;
- double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
+ double dm = mdot * dtOverride;
+ double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride;
Mass += dm;
InternalEnergy += dE;
diff --git a/Core/OrificeBoundary.cs b/Core/OrificeBoundary.cs
index 733eba2..fcf8785 100644
--- a/Core/OrificeBoundary.cs
+++ b/Core/OrificeBoundary.cs
@@ -1,5 +1,5 @@
using System;
-using FluidSim.Components;
+using FluidSim.Interfaces;
namespace FluidSim.Core
{
diff --git a/Core/Simulation.cs b/Core/Simulation.cs
index 1e8d488..308f058 100644
--- a/Core/Simulation.cs
+++ b/Core/Simulation.cs
@@ -1,5 +1,6 @@
using System;
using FluidSim.Components;
+using FluidSim.Interfaces;
using FluidSim.Utils;
namespace FluidSim.Core
@@ -7,76 +8,64 @@ namespace FluidSim.Core
public static class Simulation
{
private static Solver solver;
- private static Volume0D volA, volB;
private static Pipe1D pipe;
- private static Connection connA, connB;
private static int stepCount;
private static double time;
private static double dt;
+ private static float sample;
+ private static double ambientPressure = 1.0 * Units.atm;
public static void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
- double V = 5.0 * Units.L;
- volA = new Volume0D(V, 2.0 * Units.atm, Units.Celsius(20), sampleRate);
- volB = new Volume0D(V, 1.0 * Units.atm, Units.Celsius(20), sampleRate);
+ double length = 0.2;
+ double radius = 5 * Units.mm;
+ double area = Units.AreaFromDiameter(radius);
- double length = 150 * Units.mm;
- double diameter = 25 * Units.mm;
- double area = Units.AreaFromDiameter(25, Units.mm);
-
- pipe = new Pipe1D(length, area, sampleRate);
- pipe.SetUniformState(volA.Density, 0.0, volA.Pressure);
- pipe.FrictionFactor = 0.02;
-
- // Connections with orifice area equal to pipe area (flange joint)
- connA = new Connection(volA.Port, pipe.PortA) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
- connB = new Connection(pipe.PortB, volB.Port) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
+ pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: 80);
+ pipe.SetUniformState(1.225, 0.0, ambientPressure);
+ pipe.FrictionFactor = 0.0;
solver = new Solver();
- solver.AddVolume(volA);
- solver.AddVolume(volB);
+ solver.SetTimeStep(dt);
solver.AddPipe(pipe);
- solver.AddConnection(connA);
- solver.AddConnection(connB);
+ solver.SetPipeBoundary(pipe, isLeft: true, BoundaryType.OpenEnd, ambientPressure);
+ solver.SetPipeBoundary(pipe, isLeft: false, BoundaryType.ClosedEnd);
+
+ // Excite the pipe with an initial pressure pulse near the open end
+ int pulseCells = 5;
+ double pulsePressure = 4 * ambientPressure;
+ for (int i = 0; i < pulseCells; i++)
+ pipe.SetCellState(i, 1.225, 0.0, pulsePressure);
}
public static float Process()
{
- solver.Step();
+ sample = solver.Step();
time += dt;
stepCount++;
+
+ // Override the audio sample with mid-pipe pressure deviation
+ double pMid = pipe.GetPressureAtFraction(0.5);
+ sample = (float)((pMid - ambientPressure) / ambientPressure);
+
Log();
- return 0f;
+ return sample;
}
public static void Log()
{
- bool logPipe = true;
-
- if ((stepCount <= 10 || (stepCount <= 1000 && stepCount % 100 == 0)) || stepCount % 1000 == 0 && stepCount < 10000)
+ if (stepCount % 10 == 0 && stepCount < 1000)
{
- // Summary line
+ double pMid = pipe.GetPressureAtFraction(0.5);
+ double pOpen = pipe.GetCellPressure(0);
+ double pClosed = pipe.GetCellPressure(pipe.GetCellCount() - 1);
Console.WriteLine(
$"t = {time * 1e3:F3} ms Step {stepCount:D4}: " +
- $"PA = {volA.Pressure / 1e5:F6} bar, " +
- $"PB = {volB.Pressure / 1e5:F6} bar, " +
- $"FlowA = {pipe.PortA.MassFlowRate * 1e3:F2} g/s");
-
- // Per‑cell state
- if (logPipe && stepCount <= 1000)
- {
- int n = pipe.GetCellCount();
- for (int i = 0; i < n; i++)
- {
- double rho = pipe.GetCellDensity(i);
- double p = pipe.GetCellPressure(i);
- double u = pipe.GetCellVelocity(i);
- Console.WriteLine(
- $" Cell {i,2}: ρ={rho,8:F4} kg/m³, p={p,10:F2} Pa, u={u,8:F3} m/s");
- }
- }
+ $"Sample: = {sample:F3}, " +
+ $"P_mid = {pMid:F2} Pa ({pMid / ambientPressure:F4} atm), " +
+ $"P_open = {pOpen:F2} Pa, P_closed = {pClosed:F2} Pa");
}
}
}
diff --git a/Core/Solver.cs b/Core/Solver.cs
index 0c7807f..5b5b2c6 100644
--- a/Core/Solver.cs
+++ b/Core/Solver.cs
@@ -10,70 +10,163 @@ namespace FluidSim.Core
private readonly List _pipes = new();
private readonly List _connections = new();
+ private double _dt; // global time step
+
public void AddVolume(Volume0D v) => _volumes.Add(v);
public void AddPipe(Pipe1D p) => _pipes.Add(p);
public void AddConnection(Connection c) => _connections.Add(c);
- public void Step()
+ /// Set the global time step (called from Simulation).
+ public void SetTimeStep(double dt) => _dt = dt;
+
+ ///
+ /// Convenient method to set the boundary type of a pipe end.
+ ///
+ public void SetPipeBoundary(Pipe1D pipe, bool isLeft, BoundaryType type, double ambientPressure = 101325.0)
{
- // 1. Volumes publish state to their ports
+ if (isLeft)
+ {
+ pipe.SetLeftBoundaryType(type);
+ if (type == BoundaryType.OpenEnd)
+ pipe.SetLeftAmbientPressure(ambientPressure);
+ }
+ else
+ {
+ pipe.SetRightBoundaryType(type);
+ if (type == BoundaryType.OpenEnd)
+ pipe.SetRightAmbientPressure(ambientPressure);
+ }
+ }
+
+ public float Step()
+ {
+ // 1. Volumes publish state to ports (only needed if any volume exists)
foreach (var v in _volumes)
v.PushStateToPort();
- // 2. Set volume states as boundary conditions on pipes
+ // 2. Set initial pipe boundary conditions ONLY for volume‑coupled ends
foreach (var conn in _connections)
{
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
- SetVolumeBC(conn.PortA, conn.PortB);
+ {
+ var pipe = GetPipe(conn.PortA);
+ bool isLeft = pipe.PortA == conn.PortA;
+ BoundaryType bc = isLeft ? pipe.LeftBCType : pipe.RightBCType;
+ if (bc == BoundaryType.VolumeCoupling)
+ SetVolumeBC(conn.PortA, conn.PortB);
+ }
else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
- SetVolumeBC(conn.PortB, conn.PortA);
+ {
+ var pipe = GetPipe(conn.PortB);
+ bool isLeft = pipe.PortB == conn.PortB;
+ BoundaryType bc = isLeft ? pipe.LeftBCType : pipe.RightBCType;
+ if (bc == BoundaryType.VolumeCoupling)
+ SetVolumeBC(conn.PortB, conn.PortA);
+ }
}
- // 3. Run pipe simulations
+ // 3. Determine number of sub‑steps
+ int nSub = 1;
foreach (var p in _pipes)
- p.Simulate();
+ nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
+ double dtSub = _dt / nSub;
- // 4. Transfer pipe‑port flows to volume ports
- foreach (var conn in _connections)
+ // 4. Sub‑step loop
+ for (int sub = 0; sub < nSub; sub++)
{
- if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
- TransferPipeToVolume(conn.PortA, conn.PortB);
- else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
- TransferPipeToVolume(conn.PortB, conn.PortA);
+ foreach (var p in _pipes)
+ p.SimulateSingleStep(dtSub);
+
+ // Transfer flows only for volume‑coupled connections
+ foreach (var conn in _connections)
+ {
+ if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+ {
+ var pipe = GetPipe(conn.PortA);
+ bool isLeft = pipe.PortA == conn.PortA;
+ if (pipe.LeftBCType == BoundaryType.VolumeCoupling || pipe.RightBCType == BoundaryType.VolumeCoupling)
+ TransferAndIntegrate(conn.PortA, conn.PortB, dtSub);
+ }
+ else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+ {
+ var pipe = GetPipe(conn.PortB);
+ bool isLeft = pipe.PortB == conn.PortB;
+ if (pipe.LeftBCType == BoundaryType.VolumeCoupling || pipe.RightBCType == BoundaryType.VolumeCoupling)
+ TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
+ }
+ }
+
+ // Update BCs for volume‑coupled ends between sub‑steps
+ if (sub < nSub - 1)
+ {
+ foreach (var v in _volumes)
+ v.PushStateToPort();
+
+ foreach (var conn in _connections)
+ {
+ if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+ {
+ var pipe = GetPipe(conn.PortA);
+ bool isLeft = pipe.PortA == conn.PortA;
+ if ((isLeft && pipe.LeftBCType == BoundaryType.VolumeCoupling) ||
+ (!isLeft && pipe.RightBCType == BoundaryType.VolumeCoupling))
+ SetVolumeBC(conn.PortA, conn.PortB);
+ }
+ else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+ {
+ var pipe = GetPipe(conn.PortB);
+ bool isLeft = pipe.PortB == conn.PortB;
+ if ((isLeft && pipe.LeftBCType == BoundaryType.VolumeCoupling) ||
+ (!isLeft && pipe.RightBCType == BoundaryType.VolumeCoupling))
+ SetVolumeBC(conn.PortB, conn.PortA);
+ }
+ }
+ }
}
- // 5. Integrate volumes
- foreach (var v in _volumes)
- v.Integrate();
+ // 5. Audio samples from SoundConnections (if any)
+ var audioSamples = new List();
+ foreach (var conn in _connections)
+ {
+ if (conn is SoundConnection sc)
+ audioSamples.Add(sc.GetAudioSample());
+ }
+
+ // 6. Clear volume BC flags
+ foreach (var p in _pipes)
+ p.ClearBC();
+
+ return SoundProcessor.MixAndClip(audioSamples.ToArray());
}
- bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
- bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
- Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
+ private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
+ private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
+ private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
+ private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p);
- void SetVolumeBC(Port pipePort, Port volPort)
+ private void SetVolumeBC(Port pipePort, Port volPort)
{
- Pipe1D pipe = GetPipe(pipePort);
+ var pipe = GetPipe(pipePort);
if (pipe == null) return;
bool isLeft = pipe.PortA == pipePort;
-
if (isLeft)
pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
else
pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
}
- void TransferPipeToVolume(Port pipePort, Port volPort)
+ private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
{
double mdot = pipePort.MassFlowRate;
volPort.MassFlowRate = -mdot;
if (mdot < 0) // pipe → volume
{
- // pipePort.SpecificEnthalpy is already total (h + ½u²)
volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
}
- // else: volume → pipe, volume’s own static enthalpy is used (already set)
+ // else: volume’s own enthalpy (set by PushStateToPort) is used
+
+ GetVolume(volPort)?.Integrate(dtSub);
}
}
}
\ No newline at end of file
diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs
new file mode 100644
index 0000000..d67d54f
--- /dev/null
+++ b/Core/SoundProcessor.cs
@@ -0,0 +1,23 @@
+namespace FluidSim.Core
+{
+ ///
+ /// Mixes multiple audio samples and applies a soft‑clipping tanh.
+ ///
+ public static class SoundProcessor
+ {
+ /// Overall gain applied after mixing (before tanh).
+ public static float MasterGain { get; set; } = 0.01f;
+
+ ///
+ /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1].
+ ///
+ public static float MixAndClip(params float[] samples)
+ {
+ float sum = 0f;
+ foreach (float s in samples)
+ sum += s;
+ sum *= MasterGain;
+ return sum;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/Connection.cs b/Interfaces/Connection.cs
similarity index 88%
rename from Components/Connection.cs
rename to Interfaces/Connection.cs
index 34dcedd..409542c 100644
--- a/Components/Connection.cs
+++ b/Interfaces/Connection.cs
@@ -1,6 +1,4 @@
-using FluidSim.Interfaces;
-
-namespace FluidSim.Components
+namespace FluidSim.Interfaces
{
/// Pure data link between two ports, with orifice parameters.
public class Connection
diff --git a/Interfaces/SoundConnection.cs b/Interfaces/SoundConnection.cs
new file mode 100644
index 0000000..304ce7a
--- /dev/null
+++ b/Interfaces/SoundConnection.cs
@@ -0,0 +1,25 @@
+namespace FluidSim.Interfaces
+{
+ ///
+ /// A Connection that also produces an audio sample from the pressure drop across it.
+ ///
+ public class SoundConnection : Connection
+ {
+ /// Gain applied to the normalised pressure difference.
+ public float Gain { get; set; } = 1.0f;
+
+ /// Reference pressure used for normalisation (Pa). Default: 1 atm.
+ public double ReferencePressure { get; set; } = 101325.0;
+
+ public SoundConnection(Port a, Port b) : base(a, b) { }
+
+ ///
+ /// Returns a normalised audio sample proportional to the pressure difference.
+ ///
+ public float GetAudioSample()
+ {
+ double dp = PortA.Pressure - PortB.Pressure;
+ return (float)(dp / ReferencePressure) * Gain;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Utils/Units.cs b/Utils/Units.cs
index b151478..3891b1e 100644
--- a/Utils/Units.cs
+++ b/Utils/Units.cs
@@ -21,10 +21,10 @@ namespace FluidSim.Utils
public static double Celsius(double tC) => tC + 273.15;
- public static double AreaFromRadius(double radius, double unit = mm) =>
- Math.PI * (radius * unit) * (radius * unit);
+ public static double AreaFromRadius(double radius) =>
+ Math.PI * (radius) * (radius);
- public static double AreaFromDiameter(double diameter, double unit = mm) =>
- Math.PI * 0.25 * (diameter * unit) * (diameter * unit);
+ public static double AreaFromDiameter(double diameter) =>
+ Math.PI * 0.25 * (diameter) * (diameter);
}
}
\ No newline at end of file