Open end working

Components/Pipe1D.cs

@@ -4,10 +4,8 @@ using FluidSim.Interfaces;
 namespace FluidSim.Components
 {
     /// <summary>
-    /// 1‑D compressible Euler pipe (finite‑volume, HLLC flux).
+    /// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme.
-    /// Boundary conditions are set externally via SetGhostLeft/Right.
+    /// Ghost states are set externally via SetGhostLeft/Right; they are always required.
-    /// Enforces that ghosts are always valid before stepping.
-    /// Uses exponential damping and Newtonian energy relaxation.
     /// </summary>
     public class Pipe1D : IComponent
     {
@@ -17,7 +15,6 @@ namespace FluidSim.Components
         public double DampingMultiplier { get; set; } = 1.0;
         public double EnergyRelaxationRate { get; set; } = 0.0;   // 1/s
 
-        // Ambient pressure for the energy relaxation term (default 101325 Pa)
         private double _ambientPressure = 101325.0;
         public double AmbientPressure
         {
@@ -29,37 +26,21 @@ namespace FluidSim.Components
             }
         }
 
-        // Geometry
+        private readonly int _n;
-        private readonly int _n;                // number of real cells
+        private readonly double _dx;
-        private readonly double _dx;            // cell size (m)
+        private readonly double _diameter;
-        private readonly double _diameter;      // m
         private readonly double _gamma = 1.4;
 
-        // Conserved variables [0 .. _n-1]
+        private double[] _rho, _rhou, _E;
-        private double[] _rho;
+        private double[] _fluxM, _fluxP, _fluxE;   // flux at cell faces (0.._n)
-        private double[] _rhou;
-        private double[] _E;
 
-        // Face fluxes [0 .. _n]
-        private double[] _fluxM;
-        private double[] _fluxP;
-        private double[] _fluxE;
-
-        // Ghost cells (set externally)
         private double _rhoGhostL, _uGhostL, _pGhostL;
         private double _rhoGhostR, _uGhostR, _pGhostR;
         private bool _ghostLValid, _ghostRValid;
 
-        // Pre‑computed damping coefficient
         private double _laminarCoeff;
-        private double _ambientEnergyReference;   // internal energy density at ambient pressure
+        private double _ambientEnergyReference;
 
-        /// <summary>
-        /// Initialise a pipe with a given cell count.
-        /// </summary>
-        /// <param name="length">Pipe length (m).</param>
-        /// <param name="area">Cross‑sectional area (m²).</param>
-        /// <param name="cellCount">Number of finite‑volume cells (≥ 4).</param>
         public Pipe1D(double length, double area, int cellCount)
         {
             if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4");
@@ -77,48 +58,32 @@ namespace FluidSim.Components
             _fluxP = new double[_n + 1];
             _fluxE = new double[_n + 1];
 
-            // Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step)
             double mu_air = 1.8e-5;
             double radius = _diameter * 0.5;
             _laminarCoeff = 8.0 * mu_air / (radius * radius);
 
-            // Ambient energy reference (internal energy per unit volume at 101325 Pa)
             _ambientEnergyReference = 101325.0 / (_gamma - 1.0);
 
             PortA = new Port { Owner = this };
             PortB = new Port { Owner = this };
 
-            // Initial state = still air at ambient conditions
             SetUniformState(1.225, 0.0, 101325.0);
         }
 
         IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };
 
-        // No integration needed for pipes – state is advanced via sub‑steps
         public void UpdateState(double dt) { }
 
-        // ---------- Ghost cell interface ----------
+        // ---------- Ghost interface ----------
         public void SetGhostLeft(double rho, double u, double p)
         {
-            _rhoGhostL = rho;
+            _rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _ghostLValid = true;
-            _uGhostL = u;
-            _pGhostL = p;
-            _ghostLValid = true;
         }
-
         public void SetGhostRight(double rho, double u, double p)
         {
-            _rhoGhostR = rho;
+            _rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _ghostRValid = true;
-            _uGhostR = u;
-            _pGhostR = p;
-            _ghostRValid = true;
-        }
-
-        public void ClearGhostFlags()
-        {
-            _ghostLValid = false;
-            _ghostRValid = false;
         }
+        public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; }
 
         public (double rho, double u, double p) GetInteriorStateLeft()
         {
@@ -127,7 +92,6 @@ namespace FluidSim.Components
             double p = PressureScalar(0);
             return (rho, u, p);
         }
-
         public (double rho, double u, double p) GetInteriorStateRight()
         {
             double rho = Math.Max(_rho[_n - 1], 1e-12);
@@ -135,18 +99,11 @@ namespace FluidSim.Components
             double p = PressureScalar(_n - 1);
             return (rho, u, p);
         }
-
         public int CellCount => _n;
-
         public double GetCellDensity(int i) => _rho[i];
-        public double GetCellVelocity(int i)
+        public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
-        {
-            double rho = Math.Max(_rho[i], 1e-12);
-            return _rhou[i] / rho;
-        }
         public double GetCellPressure(int i) => PressureScalar(i);
 
-        // ---------- Sub‑stepping ----------
         public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
         {
             double maxW = 0.0;
@@ -163,38 +120,65 @@ namespace FluidSim.Components
             return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
         }
 
-        // ---------- Main simulation step (per sub‑step) ----------
+        // ---------- Main step (per sub‑step) ----------
         public void SimulateSingleStep(double dtSub)
         {
-            // Enforce that both ends have been provided with ghost states
             if (!_ghostLValid || !_ghostRValid)
-                throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep.");
+                throw new InvalidOperationException("Ghost cells not set before SimulateSingleStep.");
 
             double dt = dtSub;
             int n = _n;
 
-            // Left boundary face (index 0)
+            // ---- Compute fluxes at all faces using Lax‑Friedrichs ----
-            HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0),
+            // Left face (0): between ghostL and cell 0
-                     out _fluxM[0], out _fluxP[0], out _fluxE[0]);
+            double rL = Math.Max(_rhoGhostL, 1e-12);
+            double pL = _pGhostL;
+            double uL = _uGhostL;
+            double eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
+
+            double rR = Math.Max(_rho[0], 1e-12);
+            double pR = PressureScalar(0);
+            double uR = _rhou[0] / rR;
+            double eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
+
+            LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
+                              out _fluxM[0], out _fluxP[0], out _fluxE[0]);
 
             // Internal faces (1 .. n-1)
             for (int f = 1; f < n; f++)
             {
-                double rhoL = Math.Max(_rho[f - 1], 1e-12);
+                int iL = f - 1;
-                double uL = _rhou[f - 1] / rhoL;
+                int iR = f;
-                double pL = PressureScalar(f - 1);
+
-                double rhoR = Math.Max(_rho[f], 1e-12);
+                rL = Math.Max(_rho[iL], 1e-12);
-                double uR = _rhou[f] / rhoR;
+                pL = PressureScalar(iL);
-                double pR = PressureScalar(f);
+                uL = _rhou[iL] / rL;
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]);
+                eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
+
+                rR = Math.Max(_rho[iR], 1e-12);
+                pR = PressureScalar(iR);
+                uR = _rhou[iR] / rR;
+                eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
+
+                LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
+                                  out _fluxM[f], out _fluxP[f], out _fluxE[f]);
             }
 
-            // Right boundary face (index n)
+            // Right face (n): between cell n-1 and ghostR
-            HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1),
+            rL = Math.Max(_rho[n - 1], 1e-12);
-                     _rhoGhostR, _uGhostR, _pGhostR,
+            pL = PressureScalar(n - 1);
-                     out _fluxM[n], out _fluxP[n], out _fluxE[n]);
+            uL = _rhou[n - 1] / rL;
+            eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
 
-            // Cell update
+            rR = Math.Max(_rhoGhostR, 1e-12);
+            pR = _pGhostR;
+            uR = _uGhostR;
+            eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
+
+            LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
+                              out _fluxM[n], out _fluxP[n], out _fluxE[n]);
+
+            // ---- Cell update ----
             double dt_dx = dt / _dx;
             double coeff = _laminarCoeff * DampingMultiplier;
             double relaxRate = EnergyRelaxationRate;
@@ -213,15 +197,12 @@ namespace FluidSim.Components
                 double newRu = ru - dt_dx * dP;
                 double newE = E - dt_dx * dE_flux;
 
-                // Wall friction damping (laminar)
                 double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
                 newRu *= dampingFactor;
 
-                // Newtonian cooling toward ambient energy
                 double relaxFactor = Math.Exp(-relaxRate * dt);
                 newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
 
-                // Clamps – minimum density 1e-12, minimum pressure 100 Pa
                 newR = Math.Max(newR, 1e-12);
                 double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
                 double eMin = 100.0 / (_gamma - 1.0) + kin;
@@ -232,77 +213,54 @@ namespace FluidSim.Components
                 _E[i] = newE;
             }
 
-            // Update port states to reflect the current interior state (for audio / monitoring)
+            // Update port states
             (double rhoA, double uA, double pA) = GetInteriorStateLeft();
-            PortA.Pressure = pA;
+            PortA.Pressure = pA; PortA.Density = rhoA;
-            PortA.Density = rhoA;
             PortA.Temperature = pA / (rhoA * 287.0);
             PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA;
 
             (double rhoB, double uB, double pB) = GetInteriorStateRight();
-            PortB.Pressure = pB;
+            PortB.Pressure = pB; PortB.Density = rhoB;
-            PortB.Density = rhoB;
             PortB.Temperature = pB / (rhoB * 287.0);
             PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB;
         }
 
-        // ---------- Private helpers ----------
+        // ---------- Lax‑Friedrichs flux ----------
+        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)
+        {
+            // Primitive states
+            double rhoL = rL, rhoR = rR;
+            double EL = rhoL * eL;   // total energy per volume = rho * (specific total energy)
+            double ER = rhoR * eR;
+
+            // Conserved vectors U = (ρ, ρu, E)
+            // Flux F = (ρu, ρu²+p, (E+p)u)
+            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;
+
+            // Lax‑Friedrichs dissipation coefficient α = max(|u|+c) over whole domain, but here we use local max to be simple:
+            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);
             return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
         }
 
-        /// <summary>
-        /// HLLC approximate Riemann solver (Toro, 1997).
-        /// Computes the numerical flux at a face given left and right states.
-        /// </summary>
-        private 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 / rL);
-            double cR = Math.Sqrt(_gamma * pR / rR);
-            double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;   // specific total energy
-            double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
-
-            // Wave speed estimates (Davis, 1988)
-            double SL = Math.Min(uL - cL, uR - cR);
-            double SR = Math.Max(uL + cL, uR + cR);
-
-            double denom = rL * (SL - uL) - rR * (SR - uR);
-            double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
-
-            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;
-
-            if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; }
-            else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; }
-            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 = pR + rR * (SR - uR) * (Ss - uR);
-                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
-                fm = rsR * Ss;
-                fp = rsR * Ss * Ss + ps;
-                fe = (rsR * EsR + ps) * Ss;
-            }
-        }
-
-        /// <summary>Initialise all cells to a uniform state (rho, u, p).</summary>
         public void SetUniformState(double rho, double u, double p)
         {
             double e = p / ((_gamma - 1.0) * rho);
@@ -314,5 +272,24 @@ namespace FluidSim.Components
                 _E[i] = E;
             }
         }
+
+        public void SetCellState(int i, double rho, double u, double p)
+        {
+            if (i < 0 || i >= _n) return;
+            double e = p / ((_gamma - 1.0) * rho);
+            double E = rho * e + 0.5 * rho * u * u;
+            _rho[i] = rho;
+            _rhou[i] = rho * u;
+            _E[i] = E;
+        }
+
+        public void SetCellPressure(int i, double p)
+        {
+            if (i < 0 || i >= _n) return;
+            double rho = _rho[i];
+            double u = _rhou[i] / rho;
+            double e = p / ((_gamma - 1.0) * rho);
+            _E[i] = rho * e + 0.5 * rho * u * u;
+        }
     }
 }

Core/IsentropicOrifice.cs

@@ -4,38 +4,38 @@ namespace FluidSim.Core
 {
     /// <summary>
     /// Compressible flow through an orifice, modelled as an isentropic nozzle.
-    /// Supports choked and unchoked flow, forward and reverse.
+    /// The caller provides the upstream stagnation state (pUp, rhoUp, TUp),
+    /// downstream pressure, orifice area, discharge coefficient, and gas properties.
+    /// Returns the face state and mass flow from upstream to downstream.
     /// </summary>
     public static class IsentropicOrifice
     {
-        /// <summary>
+        public static void Compute(
-        /// Compute mass flow and face primitive state for an orifice.
+            double pUp, double rhoUp, double TUp,  // upstream stagnation
-        /// </summary>
+            double pDown,                           // downstream back pressure
-        /// <param name="pUp">Upstream stagnation pressure (Pa).</param>
+            double gamma, double R, double area, double Cd,
-        /// <param name="rhoUp">Upstream stagnation density (kg/m³).</param>
+            out double mdot, out double rhoFace, out double uFace, out double pFace)
-        /// <param name="gamma">Ratio of specific heats.</param>
-        /// <param name="R">Specific gas constant (J/kg·K).</param>
-        /// <param name="pDown">Downstream static pressure (Pa).</param>
-        /// <param name="area">Effective orifice area (m²).</param>
-        /// <param name="Cd">Discharge coefficient (default 0.62).</param>
-        /// <param name="mdot">Mass flow rate (kg/s), positive from upstream to downstream.</param>
-        /// <param name="rhoFace">Face density (kg/m³).</param>
-        /// <param name="uFace">Face velocity (m/s).</param>
-        /// <param name="pFace">Face pressure (Pa).</param>
-        public static void Compute(double pUp, double rhoUp, double TUp, double gamma, double R,
-                           double pDown, double area, double Cd,
-                           out double mdot, out double rhoFace, out double uFace, out double pFace)
         {
-            // mdot is positive from upstream to downstream.
+            mdot = 0; rhoFace = rhoUp; uFace = 0; pFace = pUp;
-            double pr = Math.Max(pDown / pUp, 1e-6);
-            double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
-            if (pr < prCrit) pr = prCrit;
 
-            double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
+            if (area <= 0 || pUp <= 0 || rhoUp <= 0 || TUp <= 0)
-            uFace = M * Math.Sqrt(gamma * R * TUp);
+                return;
+
+            double pr = pDown / pUp;
+            if (pr < 1e-6) pr = 1e-6;
+
+            double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
+            if (pr < prCrit) pr = prCrit;   // choked flow
+
+            double exponent = (gamma - 1.0) / gamma;
+            double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -exponent) - 1.0));
+            if (double.IsNaN(M)) M = 0;
+
+            double aUp = Math.Sqrt(gamma * R * TUp);
+            uFace = M * aUp;
             rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma);
             pFace = pUp * pr;
-            mdot = rhoFace * uFace * area * Cd;   // mass flow from upstream to downstream
+            mdot = rhoFace * uFace * area * Cd;    // positive from upstream to downstream
         }
     }
 }

Core/OpenEndLink.cs

@@ -3,19 +3,15 @@ using FluidSim.Components;
 
 namespace FluidSim.Core
 {
-    /// <summary>
-    /// Characteristic open‑end boundary condition.
-    /// For subsonic outflow the outgoing Riemann invariant is conserved,
-    /// and the ghost pressure is set to the prescribed ambient value.
-    /// </summary>
     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;
 
-        // Last resolved state (for audio / monitoring)
         public double LastMassFlowRate { get; private set; }
         public double LastFaceDensity { get; private set; }
         public double LastFaceVelocity { get; private set; }
@@ -27,9 +23,6 @@ namespace FluidSim.Core
             IsLeftEnd = isLeftEnd;
         }
 
-        /// <summary>
-        /// Compute the ghost state and mass flow for one sub‑step.
-        /// </summary>
         public void Resolve(double dtSub)
         {
             (double rhoInt, double uInt, double pInt) = IsLeftEnd
@@ -40,80 +33,61 @@ 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);
 
             double rhoGhost, uGhost, pGhost;
-            double mdot;
 
-            if (IsLeftEnd)
+            // ----- Supersonic outflow: extrapolate interior -----
+            bool supersonicOut = IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
+            if (supersonicOut)
             {
-                // Left end: outgoing invariant is J- = u - 2c/(γ-1)
+                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;
 
-                if (uInt <= -cInt) // supersonic inflow (all info from outside)
+                // Trial subsonic outflow ghost state
-                {
+                double s = pInt / Math.Pow(rhoInt, gamma);
-                    // Simple reservoir model – use ambient density and temperature 300 K
+                double rhoOut = Math.Pow(pAmb / s, 1.0 / gamma);
-                    rhoGhost = pAmb / (287.0 * 300.0);
+                double cOut = Math.Sqrt(gamma * pAmb / rhoOut);
-                    uGhost = uInt;   // keep interior velocity (should be supersonic inward)
+                double uOut = IsLeftEnd
-                    pGhost = pAmb;
+                    ? (J_minus + 2.0 * cOut / gm1)
-                }
+                    : (J_plus  - 2.0 * cOut / gm1);
-                else if (uInt < 0) // subsonic inflow
-                {
-                    double rhoAmb = pAmb / (287.0 * 300.0);
-                    double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb);
-                    uGhost = J_minus + 2.0 * cAmb / gm1;
-                    rhoGhost = rhoAmb;
-                    pGhost = pAmb;
-                }
-                else // subsonic outflow (uInt >= 0)
-                {
-                    double s = pInt / Math.Pow(rhoInt, gamma);
-                    rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma);
-                    double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost);
-                    uGhost = J_minus + 2.0 * cGhost / gm1;
-                    if (uGhost < 0) uGhost = 0;
-                    pGhost = pAmb;
-                }
-            }
-            else // Right end
-            {
-                // Right end: outgoing invariant is J+ = u + 2c/(γ-1)
-                double J_plus = uInt + 2.0 * cInt / gm1;
 
-                if (uInt >= cInt) // supersonic outflow
+                bool outflowPossible = IsLeftEnd ? (uOut <= 0) : (uOut >= 0);
+
+                if (outflowPossible)
                 {
-                    rhoGhost = rhoInt;
+                    // Subsonic outflow
-                    uGhost = uInt;
-                    pGhost = pInt;
-                }
-                else if (uInt >= 0) // subsonic outflow
-                {
-                    double s = pInt / Math.Pow(rhoInt, gamma);
-                    rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma);
-                    double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost);
-                    uGhost = J_plus - 2.0 * cGhost / gm1;
-                    if (uGhost < 0) uGhost = 0;
                     pGhost = pAmb;
+                    rhoGhost = rhoOut;
+                    uGhost = uOut;
                 }
-                else // subsonic inflow (uInt < 0)
+                else
                 {
-                    double rhoAmb = pAmb / (287.0 * 300.0);
+                    // Subsonic inflow (ambient reservoir model)
-                    double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb);
+                    pGhost = pAmb;
-                    uGhost = J_plus - 2.0 * cAmb / gm1;
                     rhoGhost = rhoAmb;
-                    pGhost = pAmb;
+                    uGhost = IsLeftEnd
+                        ? (J_minus + 2.0 * aAmb / gm1)
+                        : (J_plus  - 2.0 * aAmb / gm1);
                 }
             }
 
-            // Apply ghost to pipe
             if (IsLeftEnd)
                 Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost);
             else
                 Pipe.SetGhostRight(rhoGhost, uGhost, pGhost);
 
-            // Mass flow (positive = out of pipe)
             double area = Pipe.Area;
-            mdot = rhoGhost * uGhost * area;
+            double mdot = rhoGhost * uGhost * area;
-            if (IsLeftEnd) mdot = -mdot;   // positive u into pipe, so out of pipe is negative u
+            if (IsLeftEnd) mdot = -mdot;
             LastMassFlowRate = mdot;
             LastFaceDensity = rhoGhost;
             LastFaceVelocity = uGhost;

Core/OrificeLink.cs

@@ -6,7 +6,8 @@ namespace FluidSim.Core
 {
     /// <summary>
     /// Connects a port (volume or atmosphere) to one end of a pipe via an orifice.
-    /// The area can be dynamic (Func<double>).
+    /// Uses the isentropic nozzle model for the steady‑state relationship,
+    /// and includes acoustic inertance for dynamic (Helmholtz) behaviour.
     /// </summary>
     public class OrificeLink
     {
@@ -15,105 +16,131 @@ namespace FluidSim.Core
         public bool IsPipeLeftEnd { get; }
         public Func<double> AreaProvider { get; set; }
         public double DischargeCoefficient { get; set; } = 0.62;
-        public double Gamma { get; set; } = 1.4;
-        public double GasConstant { get; set; } = 287.0;
 
-        // Last resolved state (for audio/monitoring)
+        // Acoustic length (wall thickness + end correction) – controls the resonance frequency
+        public double EffectiveLength { get; set; } = 0.001;   // 1 mm
+
+        // Whether to include inertance; if false, uses the steady‑state nozzle model directly
+        public bool UseInertance { get; set; } = true;
+
+        // Current mass flow (kg/s, positive = volume → pipe)
+        private double _mdot;
+
         public double LastMassFlowRate { get; private set; }
         public double LastFaceDensity { get; private set; }
         public double LastFaceVelocity { get; private set; }
         public double LastFacePressure { get; private set; }
 
-        public OrificeLink(Port volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func<double> areaProvider)
+        public OrificeLink(Port? volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func<double> areaProvider)
         {
-            VolumePort = volumePort ?? throw new ArgumentNullException(nameof(volumePort));
+            VolumePort = volumePort;   // null is allowed
             Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe));
             IsPipeLeftEnd = isPipeLeftEnd;
             AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider));
+            _mdot = 0.0;
         }
 
-        /// <summary>
-        /// Resolve the coupling for one sub‑step. Computes nozzle flow (isentropic)
-        /// and sets the pipe ghost cell and the port flow rates.
-        /// </summary>
         public void Resolve(double dtSub)
         {
             double area = AreaProvider();
-            if (area < 1e-12)
+            // Closed wall or missing volume port => reflective boundary
+            if (area < 1e-12 || VolumePort == null)
             {
                 SetClosedWall();
                 return;
             }
 
-            // Retrieve volume state
+            // Gather volume state
-            double volP = VolumePort.Pressure;
+            double volP   = VolumePort.Pressure;
             double volRho = VolumePort.Density;
-            double volT = VolumePort.Temperature;
+            double volT   = VolumePort.Temperature;
-            double volH = VolumePort.SpecificEnthalpy;
+            double volH   = VolumePort.SpecificEnthalpy;
 
-            // Retrieve pipe interior state at the connected end
+            // Gather pipe interior state at the connected end
             (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd
                 ? Pipe.GetInteriorStateLeft()
                 : Pipe.GetInteriorStateRight();
 
-            // Determine upstream/downstream: if volume pressure > pipe pressure, flow is out of volume (negative into volume).
+            double pipeT = pipeP / Math.Max(pipeRho * 287.0, 1e-12);
-            bool flowOutOfVolume = volP > pipeP;
+            double gamma = 1.4;
-            double pUp, rhoUp, TUp, pDown;
+            double R     = 287.0;
-            if (flowOutOfVolume)
+
+            // ---- Steady‑state mass flow from isentropic nozzle ----
+            double mdotSS;   // positive = volume → pipe
+            double rhoFace, uFace, pFace;
+
+            if (volP >= pipeP)
             {
-                pUp = volP; rhoUp = volRho; TUp = volT; pDown = pipeP;
+                IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, DischargeCoefficient,
+                    out double mdotUpToDown, out rhoFace, out uFace, out pFace);
+                mdotSS = mdotUpToDown;   // volume → pipe
             }
             else
             {
-                // Pipe is upstream
+                IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, DischargeCoefficient,
-                pUp = pipeP; rhoUp = pipeRho; TUp = pipeP / (pipeRho * GasConstant); // temperature from pipe
+                    out double mdotUpToDown, out rhoFace, out uFace, out pFace);
-                pDown = volP;
+                mdotSS = -mdotUpToDown;  // pipe → volume → negative for volume→pipe convention
             }
 
-            // Compute isentropic nozzle flow
+            // ---- Inertance ODE (optional) ----
-            IsentropicOrifice.Compute(pUp, rhoUp, TUp, Gamma, GasConstant, pDown, area, DischargeCoefficient,
+            if (UseInertance)
-                                      out double mdotUpstreamToDown, out double rhoFace, out double uFace, out double pFace);
+            {
-
+                double rhoUp = _mdot >= 0 ? volRho : pipeRho;
-            // mdotUpstreamToDown is positive from upstream to downstream.
+                double inertance = rhoUp * EffectiveLength / area;
-            // Convert to mass flow into volume (positive mdot = into volume).
+                double dp = volP - pipeP;
-            double mdotVolume;
+                double resistance = Math.Abs(dp) / Math.Max(Math.Abs(mdotSS), 1e-12);
-            if (flowOutOfVolume)
+                double dmdot_dt = (dp - resistance * _mdot) / inertance;
-                mdotVolume = -mdotUpstreamToDown;   // out of volume is negative
+                _mdot += dmdot_dt * dtSub;
+            }
             else
-                mdotVolume = mdotUpstreamToDown;    // into volume is positive
+            {
+                _mdot = mdotSS;
+            }
 
-            // Clamp mass flow to available mass in volume (if it is a Volume0D)
+            // Clamp outflow to available mass (if finite volume)
             if (VolumePort.Owner is Volume0D vol)
             {
-                double maxMdot = vol.Mass / dtSub;
+                double maxOut = vol.Mass / dtSub;
-                if (mdotVolume > maxMdot) mdotVolume = maxMdot;
+                if (_mdot > maxOut) _mdot = maxOut;
-                if (mdotVolume < -maxMdot) mdotVolume = -maxMdot;
             }
 
-            // Apply ghost state to pipe
+            // ---- Ghost state ----
+            // Density = upstream density (consistent with current flow direction)
+            rhoFace = _mdot >= 0 ? volRho : pipeRho;
+            // Pressure = downstream pressure (consistent with nozzle exit)
+            pFace   = _mdot >= 0 ? pipeP : volP;
+            // Velocity magnitude derived from actual mass flow
+            double mdotMag = Math.Abs(_mdot);
+            uFace = mdotMag / (rhoFace * area);
+
+            if (IsPipeLeftEnd)
+                uFace = _mdot >= 0 ? uFace : -uFace;    // left end: positive u = into pipe
+            else
+                uFace = _mdot >= 0 ? -uFace : uFace;    // right end: positive u = out of pipe
+
+            // Apply ghost to pipe
             if (IsPipeLeftEnd)
                 Pipe.SetGhostLeft(rhoFace, uFace, pFace);
             else
                 Pipe.SetGhostRight(rhoFace, uFace, pFace);
 
-            // Store results
+            // ---- Store results ----
-            LastMassFlowRate = mdotVolume;
+            double mdotIntoVolume = -_mdot;   // positive = into volume
-            LastFaceDensity = rhoFace;
+            LastMassFlowRate = mdotIntoVolume;
+            LastFaceDensity  = rhoFace;
             LastFaceVelocity = uFace;
             LastFacePressure = pFace;
 
-            // Set port flow rates for volume integration
+            VolumePort.MassFlowRate = mdotIntoVolume;
-            VolumePort.MassFlowRate = mdotVolume;
+
-            if (mdotVolume >= 0)
+            // Enthalpy for volume integration
+            if (mdotIntoVolume >= 0)   // inflow → pipe enthalpy
             {
-                // Inflow: enthalpy comes from upstream (pipe)
+                double hPipe = gamma / (gamma - 1.0) * pipeP / Math.Max(pipeRho, 1e-12);
-                double pPipe = pipeP;
+                VolumePort.SpecificEnthalpy = hPipe;
-                double rhoPipe = pipeRho;
-                VolumePort.SpecificEnthalpy = Gamma / (Gamma - 1.0) * pPipe / rhoPipe;
             }
-            else
+            else                       // outflow → volume's own enthalpy
             {
-                // Outflow: volume's own specific enthalpy
                 VolumePort.SpecificEnthalpy = volH;
             }
         }
@@ -130,11 +157,12 @@ namespace FluidSim.Core
                 Pipe.SetGhostRight(rInt, -uInt, pInt);
 
             LastMassFlowRate = 0.0;
-            LastFaceDensity = rInt;
+            LastFaceDensity  = rInt;
             LastFaceVelocity = 0.0;
             LastFacePressure = pInt;
-            VolumePort.MassFlowRate = 0.0;
+            // Don't touch VolumePort if it's null
-            // Keep specific enthalpy as is (not used)
+            if (VolumePort != null)
+                VolumePort.MassFlowRate = 0.0;
         }
     }
 }

Core/SoundProcessor.cs

@@ -1,18 +1,20 @@
 using System;
-using FluidSim.Interfaces;
+using FluidSim.Core;
 
 namespace FluidSim.Core
 {
     public class SoundProcessor
     {
         private readonly double dt;
-        private readonly double scaleFactor;  // 1 / (4π r) and a user gain
+        private readonly double scaleFactor;  // 1 / (4π r)
         private double prevMassFlowOut;
-
-        // Simple low‑pass for derivative smoothing (≈ 2‑3 ms)
         private double smoothDMdt;
         private readonly double alpha;
 
+        // New: low‑pass the mass flow signal before derivative
+        private double flowLP;
+        private readonly double lpAlpha;
+
         public float Gain { get; set; } = 1.0f;
 
         public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0)
@@ -20,29 +22,34 @@ namespace FluidSim.Core
             dt = 1.0 / sampleRate;
             scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters);
 
-            // Smoothing time constant ~ 2 ms, blocks single‑sample spikes
+            // Smoothing time constant for the derivative: 10 ms (much smoother)
-            double tau = 0.002;
+            double tau = 0.010;   // 10 ms
             alpha = Math.Exp(-dt / tau);
+
+            // Low‑pass time constant for the mass flow: 5 ms (kneecap high‑freq directly)
+            double tauLP = 0.005;
+            lpAlpha = Math.Exp(-dt / tauLP);
         }
 
-        public float Process(Port port)
+        public float Process(OpenEndLink openEnd)
         {
-            // Outflow mass flow (positive = leaving pipe)
+            double flowOut = openEnd.LastMassFlowRate;
-            double flowOut = -port.MassFlowRate;
 
-            // Derivative
+            // Low‑pass the mass flow signal
-            double rawDerivative = (flowOut - prevMassFlowOut) / dt;
+            flowLP = lpAlpha * flowLP + (1.0 - lpAlpha) * flowOut;
-            prevMassFlowOut = flowOut;
 
-            // Smooth the derivative to kill isolated spikes
+            // Derivative of the smoothed mass flow
+            double rawDerivative = (flowLP - prevMassFlowOut) / dt;
+            prevMassFlowOut = flowLP;
+
+            // Smooth the derivative
             smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative;
 
             // Far‑field monopole pressure
             double pressure = smoothDMdt * scaleFactor * Gain;
 
-            // Soft clip to ±1 for audio output (safe limit)
+            // Soft clip to ±1 (should rarely trigger now)
-            float sample = (float)Math.Tanh(pressure);
+            return (float)Math.Tanh(pressure);
-            return sample;
         }
     }
 }

Program.cs

@@ -14,7 +14,8 @@ public class Program
     private static Scenario scenario;
 
     // Speed control (existing + new throttle)
-    private static double desiredSpeed = 0.01;
+    private static double desiredSpeed = 0.001;
+    //private static double desiredSpeed = 1.0;
     private static double currentSpeed = desiredSpeed;
     private const double MinSpeed = 0.0001;
     private const double MaxSpeed = 1.0;

Scenarios/TestScenario.cs

@@ -9,40 +9,41 @@ namespace FluidSim.Tests
     public class TestScenario : Scenario
     {
         private Solver solver;
-        private Volume0D volume;
         private Pipe1D pipe;
-        private Atmosphere atmosphere;
+        private OrificeLink closedEnd;    // left end – closed wall
-        private OrificeLink orificeLink;
+        private OpenEndLink openEndLink;  // right end – atmosphere
-        private OpenEndLink openEndLink;
+        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;
+        private double dt;
 
         public override void Initialize(int sampleRate)
         {
-            double dt = 1.0 / sampleRate;
+            dt = 1.0 / sampleRate;
+
+            soundProcessor = new SoundProcessor(sampleRate, 1);
+            soundProcessor.Gain = 10;
+
+            reverb = new OutdoorExhaustReverb(sampleRate);
 
             solver = new Solver();
             solver.SetTimeStep(dt);
 
-            volume = new Volume0D(1e-3, 150000.0, 300.0);
+            // Pipe: 2 m long, 1 cm² area, 200 cells
-            solver.AddComponent(volume);
+            pipe = new Pipe1D(length: 2, area: 1e-4, cellCount: 100);
-
-            pipe = new Pipe1D(2.0, 1e-4, 20);
             solver.AddComponent(pipe);
 
-            atmosphere = new Atmosphere();
+            // Initially pipe at ambient conditions
-            solver.AddComponent(atmosphere);
+            pipe.SetUniformState(1.225, 0.0, 101325.0);
 
-            // Volume → left pipe end (orifice)
+            // Left end: closed wall (area = 0 → reflective)
-            var volPort = volume.CreatePort();
+            closedEnd = new OrificeLink(null, pipe, isPipeLeftEnd: true, areaProvider: () => 0.0);
-            orificeLink = new OrificeLink(volPort, pipe, isPipeLeftEnd: true, areaProvider: () => 1e-5)
+            solver.AddOrificeLink(closedEnd);
-            {
-                DischargeCoefficient = 0.62,
-                Gamma = volume.Gamma,
-                GasConstant = volume.GasConstant
-            };
-            solver.AddOrificeLink(orificeLink);
 
-            // Right pipe end → atmosphere (characteristic open‑end)
+            // Right end: open to atmosphere
             openEndLink = new OpenEndLink(pipe, isLeftEnd: false)
             {
                 AmbientPressure = 101325.0,
@@ -51,45 +52,82 @@ namespace FluidSim.Tests
             solver.AddOpenEndLink(openEndLink);
 
             stepCount = 0;
-            Console.WriteLine("TestScenario initialized with sampleRate = " + sampleRate);
+            simTime = 0.0;
+            nextPulseTime = pulseInterval;   // first pulse at 100 ms
+
+            Console.WriteLine("Pulse reflection test – closed left, open right");
+            Console.WriteLine("Pulse injected every 100 ms at left end (cell 0)");
         }
 
         public override float Process()
         {
             solver.Step();
             stepCount++;
+            simTime += dt;
 
-            if (stepCount % 100 == 0)
+            // ---- Inject a pressure pulse at the closed end every 100 ms ----
+            if (simTime >= nextPulseTime)
             {
-                double volPressure = volume.Pressure;
+                // Apply a Gaussian pulse to the first few cells
-                double volMass = volume.Mass;
+                double ambientPressure = 101325.0;
-                double pipeLeftPressure = pipe.GetCellPressure(0);
+                double pulseAmplitude = 20 * ambientPressure;   // 0.5 atm overpressure
-                double pipeRightPressure = pipe.GetCellPressure(pipe.CellCount - 1);
+                double pulseWidth = 0.05;   // m (spread over a few cells)
-                double mdotOrifice = orificeLink.LastMassFlowRate;
+                int n = pipe.CellCount;
-                double mdotOpen = openEndLink.LastMassFlowRate;
+                double dx = 2.0 / n;
 
-                Console.WriteLine($"Step {stepCount}:");
+                // Only modify cells within 2*pulseWidth from the left end
-                Console.WriteLine($"  Vol Pressure = {volPressure:F1} Pa, Mass = {volMass:E4} kg");
+                int maxCell = Math.Min(5, n - 1);   // at most the first 5 cells
-                Console.WriteLine($"  Pipe left P = {pipeLeftPressure:F1} Pa, right P = {pipeRightPressure:F1} Pa");
+                for (int i = 0; i <= maxCell; i++)
-                Console.WriteLine($"  Orifice mdot = {mdotOrifice:E4} kg/s, Open‑end mdot = {mdotOpen:E4} kg/s");
+                {
-                Console.WriteLine();
+                    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);
+                }
+                Console.WriteLine($"Pulse injected at t = {simTime:F3} s");
+                nextPulseTime += pulseInterval;
             }
 
-            // Audio sample from the open‑end mass flow
+            // Audio from open‑end mass flow
-            return (float)openEndLink.LastMassFlowRate;
+            float sample = soundProcessor.Process(openEndLink);
+
+            // Log every 200 steps
+            if (stepCount % 1000 == 0)
+            {
+                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);
+
+                Console.WriteLine($"Step {stepCount}: P_left={pL:F1} Pa, P_mid={pM:F1} Pa, P_right={pR:F1} Pa");
+            }
+
+            if (double.IsNaN(pipe.GetCellPressure(0)))
+            {
+                Console.WriteLine("NaN detected – stopping simulation.");
+                return 0f;
+            }
+
+            return reverb.Process(sample);
         }
 
         public override void Draw(RenderWindow target)
         {
             float winWidth = target.GetView().Size.X;
             float winHeight = target.GetView().Size.Y;
-
             float pipeCenterY = winHeight / 2f;
             float margin = 60f;
             float pipeStartX = margin;
             float pipeEndX = winWidth - margin;
-
-            // Use the shared pipe drawing from the base class
             DrawPipe(target, pipe, pipeCenterY, pipeStartX, pipeEndX);
         }
     }