diff --git a/Components/Connection.cs b/Components/Connection.cs
new file mode 100644
index 0000000..34dcedd
--- /dev/null
+++ b/Components/Connection.cs
@@ -0,0 +1,17 @@
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+ /// Pure data link between two ports, with orifice parameters.
+ public class Connection
+ {
+ public Port PortA { get; }
+ public Port PortB { get; }
+
+ public double Area { get; set; } = 1e-5; // effective orifice area (m²)
+ public double DischargeCoefficient { get; set; } = 0.62;
+ public double Gamma { get; set; } = 1.4;
+
+ public Connection(Port a, Port b) => (PortA, PortB) = (a, b);
+ }
+}
\ No newline at end of file
diff --git a/Components/Orifice.cs b/Components/Orifice.cs
new file mode 100644
index 0000000..f3639e1
--- /dev/null
+++ b/Components/Orifice.cs
@@ -0,0 +1,73 @@
+using System;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+ public class Orifice
+ {
+ public Port PortA { get; }
+ public Port PortB { get; }
+
+ public double Area { get; set; }
+ public double DischargeCoeff { get; set; } = 0.62;
+ public double Gamma { get; set; } = 1.4;
+
+ public Orifice(double area)
+ {
+ Area = area;
+ PortA = new Port();
+ PortB = new Port();
+ }
+
+ public void Simulate()
+ {
+ double pA = PortA.Pressure, pB = PortB.Pressure;
+ double dp = pA - pB;
+ double rho = dp >= 0 ? PortA.Density : PortB.Density;
+ if (rho <= 0) rho = 1.225;
+
+ double massFlow;
+ double absDp = Math.Abs(dp);
+ double critical = 1e-3 * pA; // blend threshold
+
+ if (absDp < critical)
+ {
+ // Linearised region for numerical stability
+ massFlow = Area * DischargeCoeff * Math.Sqrt(2 * rho * critical) * dp / critical;
+ }
+ else
+ {
+ double sign = Math.Sign(dp);
+ double pratio = Math.Min(pA, pB) / Math.Max(pA, pB);
+ double choked = Math.Pow(2.0 / (Gamma + 1.0), Gamma / (Gamma - 1.0));
+
+ if (pratio < choked)
+ {
+ double term = Math.Sqrt(Gamma * Math.Pow(2.0 / (Gamma + 1.0), (Gamma + 1.0) / (Gamma - 1.0)));
+ massFlow = DischargeCoeff * Area * Math.Sqrt(rho * Math.Max(pA, pB)) * term;
+ }
+ else
+ {
+ double exp = 1.0 - Math.Pow(pratio, (Gamma - 1.0) / Gamma);
+ massFlow = DischargeCoeff * Area *
+ Math.Sqrt(2.0 * rho * Math.Max(pA, pB) * (Gamma / (Gamma - 1.0)) * pratio * pratio * exp);
+ }
+ massFlow *= sign;
+ }
+
+ PortA.MassFlowRate = -massFlow; // outflow from A
+ PortB.MassFlowRate = massFlow; // inflow to B
+
+ if (massFlow > 0) // A->B
+ {
+ PortA.SpecificEnthalpy = PortA.SpecificEnthalpy;
+ PortB.SpecificEnthalpy = PortA.SpecificEnthalpy;
+ }
+ else
+ {
+ PortA.SpecificEnthalpy = PortB.SpecificEnthalpy;
+ PortB.SpecificEnthalpy = PortB.SpecificEnthalpy;
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
new file mode 100644
index 0000000..027af15
--- /dev/null
+++ b/Components/Pipe1D.cs
@@ -0,0 +1,170 @@
+using System;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+ public class Pipe1D
+ {
+ public Port PortA { get; }
+ public Port PortB { get; }
+ public double Area => _area;
+
+ private int _n;
+ private double _dx, _dt, _gamma = 1.4, _area;
+ private double[] _rho, _rhou, _E;
+
+ // Boundary fluxes (set by solver before each step)
+ private double _fxL_mass, _fxL_mom, _fxL_ener;
+ private double _fxR_mass, _fxR_mom, _fxR_ener;
+ private bool _leftSet, _rightSet;
+
+ public double FrictionFactor { get; set; } = 0.02;
+
+ public int GetCellCount() => _n;
+ public double GetCellDensity(int i) => _rho[i];
+ public double GetCellPressure(int i) => Pressure(i);
+ public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
+
+ public Pipe1D(double length, double area, int nCells, int sampleRate)
+ {
+ _n = nCells;
+ _dx = length / nCells;
+ _dt = 1.0 / sampleRate;
+ _area = area;
+
+ _rho = new double[_n];
+ _rhou = new double[_n];
+ _E = new double[_n];
+
+ PortA = new Port();
+ PortB = new Port();
+ }
+
+ public void SetUniformState(double rho, double u, double p)
+ {
+ double e = p / ((_gamma - 1) * rho);
+ double Etot = rho * e + 0.5 * rho * u * u;
+ for (int i = 0; i < _n; i++)
+ {
+ _rho[i] = rho;
+ _rhou[i] = rho * u;
+ _E[i] = Etot;
+ }
+ }
+
+ public double GetLeftPressure() => Pressure(0);
+ public double GetRightPressure() => Pressure(_n - 1);
+ public double GetLeftDensity() => _rho[0];
+ public double GetRightDensity() => _rho[_n - 1];
+
+ public void SetLeftBoundaryFlux(double m, double p, double e)
+ {
+ _fxL_mass = m; _fxL_mom = p; _fxL_ener = e; _leftSet = true;
+ }
+
+ public void SetRightBoundaryFlux(double m, double p, double e)
+ {
+ _fxR_mass = m; _fxR_mom = p; _fxR_ener = e; _rightSet = true;
+ }
+
+ public void Simulate()
+ {
+ int n = _n;
+ double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];
+
+ // Left face
+ if (_leftSet) { Fm[0] = _fxL_mass; Fp[0] = _fxL_mom; Fe[0] = _fxL_ener; }
+ else { Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; } // reflective wall
+
+ // Internal faces (HLLC)
+ 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 face
+ if (_rightSet) { Fm[n] = _fxR_mass; Fp[n] = _fxR_mom; Fe[n] = _fxR_ener; }
+ else { Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; }
+
+ // Update cells
+ 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] -= _dt * dM;
+ _rhou[i] -= _dt * dP;
+ _E[i] -= _dt * dE;
+
+ // Clamp to physical
+ if (_rho[i] < 1e-12) _rho[i] = 1e-12;
+ double u = _rhou[i] / _rho[i];
+ double kinetic = 0.5 * _rho[i] * u * u;
+ if (_E[i] < kinetic) _E[i] = kinetic;
+ }
+
+ // Friction
+ if (FrictionFactor > 0)
+ {
+ double D = 2.0 * Math.Sqrt(_area / Math.PI);
+ for (int i = 0; i < n; i++)
+ {
+ double u = _rhou[i] / Math.Max(_rho[i], 1e-12);
+ double f = FrictionFactor / (2.0 * D) * _rho[i] * u * Math.Abs(u);
+ _rhou[i] -= _dt * f;
+ if (_E[i] > _dt * f * u) _E[i] -= _dt * f * u;
+ }
+ }
+
+ // Write port flows for the solver
+ PortA.MassFlowRate = _leftSet ? _fxL_mass * _area : 0.0;
+ PortB.MassFlowRate = _rightSet ? -_fxR_mass * _area : 0.0;
+
+ // Enthalpy for upwinding
+ PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * Pressure(0) / Math.Max(_rho[0], 1e-12);
+ PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * Pressure(_n - 1) / Math.Max(_rho[_n - 1], 1e-12);
+
+ // Reset for next step
+ _leftSet = _rightSet = false;
+ }
+
+ double Pressure(int i) =>
+ (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
+
+ void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
+ out double fm, out double fp, out double fe)
+ {
+ double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12));
+ double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12));
+ double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
+ double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR;
+ double SL = Math.Min(uL - cL, uR - cR);
+ double SR = Math.Max(uL + cL, uR + cR);
+ double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
+ / (rL * (SL - uL) - rR * (SR - uR));
+
+ double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
+ double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
+
+ if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
+ else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
+ else if (Ss >= 0)
+ {
+ double rsL = rL * (SL - uL) / (SL - Ss);
+ double ps = pL + rL * (SL - uL) * (Ss - uL);
+ double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
+ fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
+ }
+ else
+ {
+ double rsR = rR * (SR - uR) / (SR - Ss);
+ double ps = pL + rL * (SL - uL) * (Ss - uL);
+ double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
+ fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs
new file mode 100644
index 0000000..b5949a2
--- /dev/null
+++ b/Components/Volume0D.cs
@@ -0,0 +1,71 @@
+using System;
+using FluidSim.Interfaces;
+using FluidSim.Utils;
+
+namespace FluidSim.Components
+{
+ public class Volume0D
+ {
+ public Port Port { get; private set; }
+
+ public double Mass { get; private set; }
+ public double InternalEnergy { get; private set; }
+
+ public double Gamma { get; set; } = 1.4;
+ public double GasConstant { get; set; } = 287.0;
+
+ public double Volume { get; set; }
+ public double dVdt { get; set; }
+
+ private double _dt;
+
+ public double Density => Mass / Volume;
+ public double Pressure => (Gamma - 1.0) * InternalEnergy / Volume;
+ public double Temperature => Pressure / (Density * GasConstant);
+ public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
+
+ public Volume0D(double initialVolume, double initialPressure,
+ double initialTemperature, int sampleRate,
+ double gasConstant = 287.0, double gamma = 1.4)
+ {
+ GasConstant = gasConstant;
+ Gamma = gamma;
+ Volume = initialVolume;
+ dVdt = 0.0;
+ _dt = 1.0 / sampleRate;
+
+ double rho0 = initialPressure / (GasConstant * initialTemperature);
+ Mass = rho0 * Volume;
+ InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
+
+ Port = new Port();
+ PushStateToPort();
+ }
+
+ public void PushStateToPort()
+ {
+ Port.Pressure = Pressure;
+ Port.Density = Density;
+ Port.Temperature = Temperature;
+ Port.SpecificEnthalpy = SpecificEnthalpy;
+ }
+
+ public void Integrate()
+ {
+ double mdot = Port.MassFlowRate;
+ double h_in = Port.SpecificEnthalpy;
+
+ double dm = mdot * _dt;
+ double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
+
+ Mass += dm;
+ InternalEnergy += dE;
+
+ // Hard physical bounds – prevent NaN and unphysical states
+ if (Mass < 1e-12) Mass = 1e-12;
+ if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
+
+ PushStateToPort();
+ }
+ }
+}
\ No newline at end of file
diff --git a/Core/OrificeBoundary.cs b/Core/OrificeBoundary.cs
new file mode 100644
index 0000000..c4c22c2
--- /dev/null
+++ b/Core/OrificeBoundary.cs
@@ -0,0 +1,95 @@
+using System;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+ public static class OrificeBoundary
+ {
+ public static double MassFlow(double pA, double rhoA, double pB, double rhoB,
+ Connection conn)
+ {
+ if (double.IsNaN(pA) || double.IsNaN(rhoA) || double.IsNaN(pB) || double.IsNaN(rhoB) ||
+ double.IsInfinity(pA) || double.IsInfinity(rhoA) || double.IsInfinity(pB) || double.IsInfinity(rhoB) ||
+ pA <= 0 || rhoA <= 0 || pB <= 0 || rhoB <= 0)
+ return 0.0;
+
+ double dp = pA - pB;
+ double sign = Math.Sign(dp);
+ double absDp = Math.Abs(dp);
+ double rhoUp = dp >= 0 ? rhoA : rhoB;
+ double pUp = dp >= 0 ? pA : pB;
+ double pDown = dp >= 0 ? pB : pA;
+ double delta = 1e-6 * pUp;
+
+ if (absDp < delta)
+ {
+ double k = conn.DischargeCoefficient * conn.Area * Math.Sqrt(2 * rhoUp / delta);
+ return k * dp;
+ }
+ else
+ {
+ double pr = pDown / pUp;
+ double choked = Math.Pow(2.0 / (conn.Gamma + 1.0), conn.Gamma / (conn.Gamma - 1.0));
+ if (pr < choked)
+ {
+ double term = Math.Sqrt(conn.Gamma *
+ Math.Pow(2.0 / (conn.Gamma + 1.0), (conn.Gamma + 1.0) / (conn.Gamma - 1.0)));
+ double flow = conn.DischargeCoefficient * conn.Area *
+ Math.Sqrt(rhoUp * pUp) * term;
+ return sign * flow;
+ }
+ else
+ {
+ double ex = 1.0 - Math.Pow(pr, (conn.Gamma - 1.0) / conn.Gamma);
+ double flow = conn.DischargeCoefficient * conn.Area *
+ Math.Sqrt(2.0 * rhoUp * pUp * (conn.Gamma / (conn.Gamma - 1.0)) *
+ pr * pr * ex);
+ return sign * flow;
+ }
+ }
+ }
+
+ public static void PipeVolumeFlux(double pPipe, double rhoPipe, double uPipe,
+ double pVol, double rhoVol, double uVol,
+ Connection conn, double pipeArea,
+ bool isLeftBoundary,
+ out double massFlux, out double momFlux, out double energyFlux)
+ {
+ // mass flow from pipe to volume (positive = pipe → volume)
+ double mdot = MassFlow(pPipe, rhoPipe, pVol, rhoVol, conn);
+
+ // Limit mass flow to the amount that can leave/enter the pipe cell
+ double maxMdot = rhoPipe * pipeArea * 343.0;
+ if (Math.Abs(mdot) > maxMdot) mdot = Math.Sign(mdot) * maxMdot;
+
+ bool flowLeavesPipe = mdot > 0;
+ double uFace, pFace, rhoFace;
+ double massFluxPerArea;
+
+ if (isLeftBoundary)
+ {
+ massFluxPerArea = -mdot / pipeArea;
+ if (flowLeavesPipe)
+ { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
+ else
+ { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
+ }
+ else // right boundary
+ {
+ massFluxPerArea = mdot / pipeArea;
+ if (flowLeavesPipe)
+ { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
+ else
+ { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
+ }
+
+ // Total enthalpy of the injected fluid (corrected: mass flux × total enthalpy)
+ double specificEnthalpy = (1.4 / (1.4 - 1.0)) * pFace / Math.Max(rhoFace, 1e-12);
+ double totalEnthalpy = specificEnthalpy + 0.5 * uFace * uFace;
+
+ massFlux = massFluxPerArea;
+ momFlux = massFluxPerArea * uFace + pFace;
+ energyFlux = massFluxPerArea * totalEnthalpy;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Core/Simulation.cs b/Core/Simulation.cs
new file mode 100644
index 0000000..45a9a3c
--- /dev/null
+++ b/Core/Simulation.cs
@@ -0,0 +1,86 @@
+using System;
+using FluidSim.Components;
+using FluidSim.Utils;
+
+namespace FluidSim.Core
+{
+ public static class Simulation
+ {
+ private static Pipe1D pipe;
+ private static Connection leftConn, rightConn; // dummy connections for orifice params
+ private static double time;
+ private static double dt;
+ private static int stepCount;
+
+ public static void Initialize(int sampleRate)
+ {
+ dt = 1.0 / sampleRate;
+
+ double length = 150 * Units.mm;
+ double diameter = 25 * Units.mm;
+ double area = Units.AreaFromDiameter(25, Units.mm);
+ int nCells = 10;
+
+ pipe = new Pipe1D(length, area, nCells, sampleRate);
+ pipe.SetUniformState(1.2, 0.0, 1.0 * Units.atm); // start at 1 atm
+ pipe.FrictionFactor = 0.02;
+
+ // Dummy connections – only used for orifice parameters
+ leftConn = new Connection(null, null) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
+ rightConn = new Connection(null, null) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
+ }
+
+ public static float Process()
+ {
+ // Fixed boundary reservoirs
+ double pLeft = 1.1 * Units.atm;
+ double rhoLeft = 1.2;
+ double uLeft = 0.0;
+
+ double pRight = 1.0 * Units.atm;
+ double rhoRight = 1.2;
+ double uRight = 0.0;
+
+ // Compute boundary fluxes via orifice model
+ OrificeBoundary.PipeVolumeFlux(
+ pipe.GetLeftPressure(), pipe.GetLeftDensity(), 0.0,
+ pLeft, rhoLeft, uLeft,
+ leftConn, pipe.Area, true,
+ out double leftMassFlux, out double leftMomFlux, out double leftEnergyFlux);
+
+ OrificeBoundary.PipeVolumeFlux(
+ pipe.GetRightPressure(), pipe.GetRightDensity(), 0.0,
+ pRight, rhoRight, uRight,
+ rightConn, pipe.Area, false,
+ out double rightMassFlux, out double rightMomFlux, out double rightEnergyFlux);
+
+ pipe.SetLeftBoundaryFlux(leftMassFlux, leftMomFlux, leftEnergyFlux);
+ pipe.SetRightBoundaryFlux(rightMassFlux, rightMomFlux, rightEnergyFlux);
+ pipe.Simulate();
+
+ time += dt;
+ stepCount++;
+ Log();
+ return 0f;
+ }
+
+ public static void Log()
+ {
+ if (stepCount <= 20 || stepCount % 50 == 0)
+ {
+ Console.WriteLine($"Step {stepCount:D4} t = {time * 1e3:F3} ms");
+ for (int i = 0; i < pipe.GetCellCount(); i++)
+ {
+ double rho = pipe.GetCellDensity(i);
+ double p = pipe.GetCellPressure(i);
+ double u = pipe.GetCellVelocity(i);
+ Console.WriteLine($" Cell {i}: ρ={rho:F4} kg/m³ p={p / 1e5:F6} bar u={u:F3} m/s");
+ }
+ double leftFlow = pipe.PortA.MassFlowRate;
+ double rightFlow = pipe.PortB.MassFlowRate;
+ Console.WriteLine($" Left flow = {leftFlow * 1e3:F4} g/s Right flow = {rightFlow * 1e3:F4} g/s");
+ Console.WriteLine();
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/Core/Solver.cs b/Core/Solver.cs
new file mode 100644
index 0000000..58c2cf6
--- /dev/null
+++ b/Core/Solver.cs
@@ -0,0 +1,94 @@
+using System.Collections.Generic;
+using FluidSim.Components;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Core
+{
+ public class Solver
+ {
+ private readonly List _volumes = new();
+ private readonly List _pipes = new();
+ private readonly List _connections = new();
+
+ 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()
+ {
+ // 1. Volumes publish state
+ foreach (var v in _volumes)
+ v.PushStateToPort();
+
+ // 2. Apply orifice boundaries to pipes
+ foreach (var conn in _connections)
+ {
+ if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+ ApplyOrifice(conn, conn.PortA, conn.PortB);
+ else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+ ApplyOrifice(conn, conn.PortB, conn.PortA);
+ else if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB))
+ VolumeToVolume(conn);
+ }
+
+ // 3. Pipes simulate
+ foreach (var p in _pipes)
+ p.Simulate();
+
+ // 4. Transfer pipe flows to connected volumes
+ foreach (var conn in _connections)
+ {
+ if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+ Transfer(conn.PortA, conn.PortB);
+ else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+ Transfer(conn.PortB, conn.PortA);
+ }
+
+ // 5. Volumes integrate
+ foreach (var v in _volumes)
+ v.Integrate();
+ }
+
+ 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);
+
+ void ApplyOrifice(Connection conn, Port pipePort, Port volPort)
+ {
+ Pipe1D pipe = GetPipe(pipePort);
+ if (pipe == null) return;
+ bool isLeft = pipe.PortA == pipePort;
+
+ double pP = isLeft ? pipe.GetLeftPressure() : pipe.GetRightPressure();
+ double rhoP = isLeft ? pipe.GetLeftDensity() : pipe.GetRightDensity();
+ double uP = 0.0;
+ double pV = volPort.Pressure, rhoV = volPort.Density, uV = 0.0;
+
+ OrificeBoundary.PipeVolumeFlux(pP, rhoP, uP, pV, rhoV, uV, conn, pipe.Area,
+ isLeft, out double mf, out double pf, out double ef);
+ if (isLeft)
+ pipe.SetLeftBoundaryFlux(mf, pf, ef);
+ else
+ pipe.SetRightBoundaryFlux(mf, pf, ef);
+ }
+
+ void VolumeToVolume(Connection conn)
+ {
+ double mdot = OrificeBoundary.MassFlow(conn.PortA.Pressure, conn.PortA.Density,
+ conn.PortB.Pressure, conn.PortB.Density, conn);
+ conn.PortA.MassFlowRate = -mdot;
+ conn.PortB.MassFlowRate = mdot;
+ if (mdot > 0)
+ conn.PortB.SpecificEnthalpy = conn.PortA.SpecificEnthalpy;
+ else if (mdot < 0)
+ conn.PortA.SpecificEnthalpy = conn.PortB.SpecificEnthalpy;
+ }
+
+ void Transfer(Port pipePort, Port volPort)
+ {
+ double mdot = pipePort.MassFlowRate;
+ volPort.MassFlowRate = -mdot;
+ volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Core/SoundEngine.cs b/Core/SoundEngine.cs
new file mode 100644
index 0000000..bab240a
--- /dev/null
+++ b/Core/SoundEngine.cs
@@ -0,0 +1,131 @@
+using SFML.Audio;
+using SFML.System;
+
+namespace FluidSim;
+
+#region Lock‑free ring buffer (unchanged)
+internal class RingBuffer
+{
+ private readonly float[] buffer;
+ private volatile int readPos;
+ private volatile int writePos;
+
+ public RingBuffer(int capacity)
+ {
+ if ((capacity & (capacity - 1)) != 0)
+ throw new ArgumentException("Capacity must be a power of two.");
+ buffer = new float[capacity];
+ }
+
+ public int Count => (writePos - readPos) & (buffer.Length - 1);
+ public int Space => (readPos - writePos - 1) & (buffer.Length - 1);
+
+ public int Write(float[] data, int count)
+ {
+ int space = Space;
+ int toWrite = Math.Min(count, space);
+ int mask = buffer.Length - 1;
+ for (int i = 0; i < toWrite; i++)
+ buffer[(writePos + i) & mask] = data[i];
+ writePos = (writePos + toWrite) & mask;
+ return toWrite;
+ }
+
+ public int Read(float[] destination, int count)
+ {
+ int available = Count;
+ int toRead = Math.Min(count, available);
+ int mask = buffer.Length - 1;
+ for (int i = 0; i < toRead; i++)
+ destination[i] = buffer[(readPos + i) & mask];
+ readPos = (readPos + toRead) & mask;
+ return toRead;
+ }
+}
+#endregion
+
+#region Stereo stream that consumes the ring buffer
+internal class RingBufferStream : SoundStream
+{
+ private readonly RingBuffer ringBuffer;
+
+ public RingBufferStream(RingBuffer buffer)
+ {
+ ringBuffer = buffer;
+ // 2 channels, 44.1 kHz, standard stereo mapping
+ Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
+ }
+
+ protected override bool OnGetData(out short[] samples)
+ {
+ const int monoBlockSize = 512; // number of mono samples we'll read
+ float[] temp = new float[monoBlockSize];
+ int read = ringBuffer.Read(temp, monoBlockSize);
+ samples = new short[monoBlockSize * 2];
+
+ if (read > 0)
+ {
+ for (int i = 0; i < read; i++)
+ {
+ float clamped = Math.Clamp(temp[i], -1f, 1f);
+ short final = (short)(clamped * short.MaxValue);
+ samples[i * 2] = final; // left
+ samples[i * 2 + 1] = final; // right
+ }
+ }
+ for (int i = read * 2; i < samples.Length; i++)
+ samples[i] = 0;
+
+ return true;
+ }
+
+ protected override void OnSeek(Time timeOffset) =>
+ throw new NotSupportedException();
+}
+#endregion
+
+#region Public sound engine API (unchanged)
+public class SoundEngine : IDisposable
+{
+ private readonly RingBuffer ringBuffer;
+ private readonly RingBufferStream stream;
+ private bool isPlaying;
+
+ public SoundEngine(int bufferCapacity = 16384)
+ {
+ ringBuffer = new RingBuffer(bufferCapacity);
+ stream = new RingBufferStream(ringBuffer);
+ }
+
+ public void Start()
+ {
+ if (isPlaying) return;
+ stream.Play();
+ isPlaying = true;
+ }
+
+ public void Stop()
+ {
+ if (!isPlaying) return;
+ stream.Stop();
+ isPlaying = false;
+ float[] drain = new float[ringBuffer.Count];
+ ringBuffer.Read(drain, drain.Length);
+ }
+
+ public int WriteSamples(float[] data, int count) =>
+ ringBuffer.Write(data, count);
+
+ public float Volume
+ {
+ get => stream.Volume;
+ set => stream.Volume = value;
+ }
+
+ public void Dispose()
+ {
+ Stop();
+ stream.Dispose();
+ }
+}
+#endregion
\ No newline at end of file
diff --git a/FluidSim.csproj b/FluidSim.csproj
new file mode 100644
index 0000000..3fda142
--- /dev/null
+++ b/FluidSim.csproj
@@ -0,0 +1,16 @@
+
+
+
+ Exe
+ net10.0
+ enable
+ enable
+ true
+ true
+
+
+
+
+
+
+
diff --git a/FluidSim.slnx b/FluidSim.slnx
new file mode 100644
index 0000000..7b36487
--- /dev/null
+++ b/FluidSim.slnx
@@ -0,0 +1,3 @@
+
+
+
diff --git a/Interfaces/Junction0.cs b/Interfaces/Junction0.cs
new file mode 100644
index 0000000..882a7d1
--- /dev/null
+++ b/Interfaces/Junction0.cs
@@ -0,0 +1,10 @@
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace FluidSim.Interfaces
+{
+ internal class Junction0
+ {
+ }
+}
diff --git a/Interfaces/Junction1.cs b/Interfaces/Junction1.cs
new file mode 100644
index 0000000..5ecc76d
--- /dev/null
+++ b/Interfaces/Junction1.cs
@@ -0,0 +1,10 @@
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace FluidSim.Interfaces
+{
+ internal class Junction1
+ {
+ }
+}
diff --git a/Interfaces/Port.cs b/Interfaces/Port.cs
new file mode 100644
index 0000000..10cedbb
--- /dev/null
+++ b/Interfaces/Port.cs
@@ -0,0 +1,20 @@
+namespace FluidSim.Interfaces
+{
+ public class Port
+ {
+ public double Pressure; // Pa
+ public double MassFlowRate; // kg/s, positive INTO the component
+ public double SpecificEnthalpy; // J/kg, enthalpy of fluid entering this port
+ public double Density; // kg/m³
+ public double Temperature; // K
+
+ public Port()
+ {
+ Pressure = 101325.0;
+ MassFlowRate = 0.0;
+ SpecificEnthalpy = 0.0;
+ Density = 1.225;
+ Temperature = 300.0;
+ }
+ }
+}
\ No newline at end of file
diff --git a/Program.cs b/Program.cs
new file mode 100644
index 0000000..938ced0
--- /dev/null
+++ b/Program.cs
@@ -0,0 +1,69 @@
+using SFML.Graphics;
+using SFML.Window;
+using SFML.System;
+using System.Diagnostics;
+using FluidSim.Core;
+
+namespace FluidSim;
+
+public class Program
+{
+ private const int SampleRate = 44100;
+ private static volatile bool running = true;
+
+ public static void Main()
+ {
+ var mode = new VideoMode(new Vector2u(1280, 720));
+ var window = new RenderWindow(mode, "Fluid Simulation");
+ window.SetVerticalSyncEnabled(true);
+ window.Closed += (_, _) => { running = false; window.Close(); };
+
+ var soundEngine = new SoundEngine(bufferCapacity: 2048);
+ soundEngine.Volume = 70;
+ soundEngine.Start();
+
+ double lastAudioTime = 0.0;
+ var stopwatch = Stopwatch.StartNew();
+
+ int warmupSamples = SampleRate / 2;
+ float[] warmup = new float[warmupSamples];
+ for (int i = 0; i < warmupSamples; i++)
+ warmup[i] = 0;
+
+ soundEngine.WriteSamples(warmup, warmupSamples);
+ lastAudioTime = stopwatch.Elapsed.TotalSeconds;
+
+ const int chunkSize = 2048;
+ float[] buffer = new float[chunkSize];
+
+ Simulation.Initialize(SampleRate);
+
+ while (window.IsOpen)
+ {
+ window.DispatchEvents();
+
+ double currentTime = stopwatch.Elapsed.TotalSeconds;
+ double elapsed = currentTime - lastAudioTime;
+ int samplesNeeded = (int)(elapsed * SampleRate);
+
+ while (samplesNeeded > 0 && running)
+ {
+ int toGenerate = Math.Min(samplesNeeded, chunkSize);
+ for (int i = 0; i < toGenerate; i++)
+ {
+ buffer[i] = Simulation.Process();
+ }
+ soundEngine.WriteSamples(buffer, toGenerate);
+ samplesNeeded -= toGenerate;
+ }
+
+ lastAudioTime = currentTime;
+
+ window.Clear(Color.Black);
+ window.Display();
+ }
+
+ soundEngine.Dispose();
+ window.Dispose();
+ }
+}
\ No newline at end of file
diff --git a/Sources/EffortSource.cs b/Sources/EffortSource.cs
new file mode 100644
index 0000000..78e186c
--- /dev/null
+++ b/Sources/EffortSource.cs
@@ -0,0 +1,10 @@
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace FluidSim.Sources
+{
+ internal class EffortSource
+ {
+ }
+}
diff --git a/Sources/FlowSource.cs b/Sources/FlowSource.cs
new file mode 100644
index 0000000..2004336
--- /dev/null
+++ b/Sources/FlowSource.cs
@@ -0,0 +1,10 @@
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace FluidSim.Sources
+{
+ internal class FlowSource
+ {
+ }
+}
diff --git a/Utils/Units.cs b/Utils/Units.cs
new file mode 100644
index 0000000..b151478
--- /dev/null
+++ b/Utils/Units.cs
@@ -0,0 +1,30 @@
+using System;
+
+namespace FluidSim.Utils
+{
+ public static class Units
+ {
+ public const double mm = 1e-3;
+ public const double cm = 1e-2;
+ public const double inch = 0.0254;
+ public const double mm2 = 1e-6;
+ public const double cm2 = 1e-4;
+ public const double mL = 1e-6;
+ public const double L = 1e-3;
+ public const double Pa = 1.0;
+ public const double kPa = 1e3;
+ public const double bar = 1e5;
+ public const double atm = 101325.0;
+ public const double psi = 6894.76;
+ public const double g = 1e-3;
+ public const double lb = 0.453592;
+
+ 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 AreaFromDiameter(double diameter, double unit = mm) =>
+ Math.PI * 0.25 * (diameter * unit) * (diameter * unit);
+ }
+}
\ No newline at end of file