working helmholtz

Core/BoundarySystem.cs

@@ -6,6 +6,12 @@ namespace FluidSim.Core
 {
     public class BoundarySystem
     {
+        // ---------- Private constants ----------
+        private const float Gamma = 1.4f;
+        private const float Gm1 = Gamma - 1f;   // 0.4
+        private const float Rgas = 287f;         // J/(kg·K)
+        private const float GammaOverGm1 = Gamma / Gm1; // 3.5
+
         public struct OrificeDesc
         {
             public Port VolumePort;
@@ -19,8 +25,8 @@ namespace FluidSim.Core
             public float EffectiveLength;
             public float CurrentMdot;       // kg/s, positive = volume → pipe
 
-            // --- Dissipative loss ---
+            // --- Loss coefficient (linear resistance) ---
-            public float LossCoefficient;   // K factor for pressure drop = K * 0.5*rho*u^2
+            public float LossCoefficient;   // N·s/m⁵  or  kg/(m⁴·s)
         }
 
         public struct OpenEndDesc
@@ -50,6 +56,7 @@ namespace FluidSim.Core
         public int OrificeCount { get; private set; }
         public int OpenEndCount { get; private set; }
 
+        // ---------- Add orifice (no inertance) ----------
         public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
                                int areaIndex, float dischargeCoeff = 1f,
                                float lossCoefficient = 0f)
@@ -69,14 +76,17 @@ namespace FluidSim.Core
             OrificeCount++;
         }
 
+        // ---------- Add orifice with inertance ----------
         public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
                                             int areaIndex, float dischargeCoeff,
                                             float effectiveLength, float lossCoefficient = 0f)
         {
+            // Reuse the base AddOrifice and then override fields
             AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient);
             ref var d = ref _orifices[OrificeCount - 1];
             d.UseInertance = true;
             d.EffectiveLength = effectiveLength;
+            d.LossCoefficient = lossCoefficient;   // store the linear resistance
         }
 
         public void AddOpenEnd(int pipeIndex, bool isLeftEnd,
@@ -112,15 +122,33 @@ namespace FluidSim.Core
             return _openEnds[openEndIndex].LastFacePressure;
         }
 
+        // ---------- Resolve all orifices ----------
         public void ResolveOrifices(float dt)
         {
             for (int i = 0; i < OrificeCount; i++)
             {
                 ref var d = ref _orifices[i];
                 float area = _orificeAreas[d.AreaIndex];
+
+                // Gather volume state
+                float volP = d.VolumePort?.Pressure ?? 101325f;
+                float volRho = d.VolumePort?.Density ?? 1.2f;
+                float volT = d.VolumePort?.Temperature ?? 300f;
+                float volH = d.VolumePort?.SpecificEnthalpy ?? 0f;
+                float volAF = d.VolumePort?.AirFraction ?? 1f;
+
+                // Gather pipe interior state
+                var (pipeRho, pipeU, pipeP) = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
+                float pipeT = pipeP / MathF.Max(pipeRho * Rgas, 1e-12f);
+                float pipeAF = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
+
+                // ---- Handle closed orifice (area ≈ 0) as a wall ----
                 if (area < 1e-12f || d.VolumePort == null)
                 {
-                    // Closed wall – reflect interior state
                     var (rInt, uInt, pInt) = d.IsLeftEnd
                         ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
                         : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
@@ -137,143 +165,126 @@ namespace FluidSim.Core
                     continue;
                 }
 
-                // Gather states
+                // ---- Preliminary isentropic solution (used for face pressure if inertance is on) ----
-                float volP = d.VolumePort.Pressure;
-                float volRho = d.VolumePort.Density;
-                float volT = d.VolumePort.Temperature;
-                float volH = d.VolumePort.SpecificEnthalpy;
-                float volAF = d.VolumePort.AirFraction;
-
-                var (pipeRho, pipeU, pipeP) = d.IsLeftEnd
-                    ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
-                    : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
-                float pipeT = pipeP / MathF.Max(pipeRho * 287f, 1e-12f);
-                float pipeAF = d.IsLeftEnd
-                    ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
-                    : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
-
-                float gamma = 1.4f, R = 287f, Cd = d.DischargeCoeff;
-
-                // --- Preliminary nozzle solution (no loss) to estimate flow direction and velocity ---
                 float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
                 if (volP >= pipeP)
                 {
-                    IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, Cd,
+                    IsentropicOrifice.Compute(volP, volRho, volT, pipeP, Gamma, Rgas, area, d.DischargeCoeff,
                         out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
                 }
                 else
                 {
-                    IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, Cd,
+                    IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, Gamma, Rgas, area, d.DischargeCoeff,
                         out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
                     mdotEst = -mdotEst;
                 }
 
-                // --- Apply symmetric loss if LossCoefficient > 0 ---
+                // ---- Compute ghost state ----
-                float volP_eff = volP;
+                float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
-                float pipeP_eff = pipeP;
-                if (d.LossCoefficient > 0f && MathF.Abs(mdotEst) > 1e-12f)
-                {
-                    float rhoRef = mdotEst >= 0 ? rhoFaceEst : rhoFaceEst; // rhoFaceEst already reflects the correct side
-                    float uRef = uFaceEst;
-                    float dynP = 0.5f * rhoRef * uRef * uRef * d.LossCoefficient;
 
-                    // Clamp the loss to avoid overshoot (max 80% of pressure difference)
-                    float dp = MathF.Abs(volP - pipeP);
-                    dynP = MathF.Min(dynP, 0.8f * dp);
-
-                    // Apply symmetrically: loss reduces the higher pressure and increases the lower pressure
-                    if (mdotEst >= 0) // volume → pipe
-                    {
-                        volP_eff -= dynP;
-                        pipeP_eff += dynP;
-                    }
-                    else               // pipe → volume
-                    {
-                        pipeP_eff -= dynP;
-                        volP_eff += dynP;
-                    }
-                }
-
-                // --- Final nozzle solution with corrected pressures ---
-                float mdotSS, rhoFace0, uFace0, pFace0;
-                if (volP_eff >= pipeP_eff)
-                {
-                    IsentropicOrifice.Compute(volP_eff, volRho, volT, pipeP_eff, gamma, R, area, Cd,
-                        out float mUp, out rhoFace0, out uFace0, out pFace0);
-                    mdotSS = mUp;
-                }
-                else
-                {
-                    IsentropicOrifice.Compute(pipeP_eff, pipeRho, pipeT, volP_eff, gamma, R, area, Cd,
-                        out float mUp, out rhoFace0, out uFace0, out pFace0);
-                    mdotSS = -mUp;
-                }
-
-                float mdot;
                 if (d.UseInertance)
                 {
+                    // ---- Inertance ODE with loss ----
                     float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
-                    float inertance = rhoUp * d.EffectiveLength / area;
+                    float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
-                    float dp = volP_eff - pipeP_eff;
+                    float dp = volP - pipeP;
-                    float resistance = MathF.Abs(dp) / MathF.Max(MathF.Abs(mdotSS), 1e-12f);
+                    float Rlin = d.LossCoefficient;   // linear resistance
-                    float dmdot_dt = (dp - resistance * d.CurrentMdot) / inertance;
-                    mdot = d.CurrentMdot + dmdot_dt * dt;
 
+                    // Forward Euler
+                    float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
+                    float mdotNew = d.CurrentMdot + dmdot_dt * dt;
+
+                    // ---------- Symmetric flow limiter ----------
+                    // 1) limit by cavity mass
                     if (d.VolumePort.Owner is Volume0D vol0)
                     {
                         float maxOut = vol0.Mass / dt;
-                        if (mdot > maxOut) mdot = maxOut;
+                        if (mdotNew > maxOut) mdotNew = maxOut;   // cavity → pipe
+                        if (mdotNew < -maxOut) mdotNew = -maxOut;   // pipe → cavity
                     }
-                    if (float.IsNaN(mdot)) mdot = 0f;
+
+                    // 2) limit by mass in the adjacent pipe cell
+                    int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
+                                              : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
+                    float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
+                    float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);   // uses the new public method
+                    float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
+                    float maxFromPipe = pipeCellMass / dt;
+                    if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;   // prevent emptying the pipe cell
+
+                    // NaN safety
+                    if (float.IsNaN(mdotNew)) mdotNew = 0f;
+
+                    // Store for next step
+                    d.CurrentMdot = mdotNew;
+                    mdotFinal = mdotNew;
+
+                    // Ghost state: use the isentropic face pressure (pFaceEst) as reference,
+                    // but compute velocity from mdotFinal.
+                    rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
+                    pFace = pFaceEst;   // approximate face pressure
+                    uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
                 }
                 else
                 {
-                    mdot = mdotSS;
+                    // ---- Standard quasi‑steady orifice ----
+                    mdotFinal = mdotEst;
+                    rhoFace = rhoFaceEst;
+                    uFace = uFaceEst;
+                    pFace = pFaceEst;
+
+                    // Limit outflow from cavity
                     if (d.VolumePort.Owner is Volume0D vol0)
                     {
                         float maxOut = vol0.Mass / dt;
-                        if (mdot > maxOut) mdot = maxOut;
+                        if (mdotFinal > maxOut) mdotFinal = maxOut;
                     }
+                    d.CurrentMdot = mdotFinal;
                 }
 
-                d.CurrentMdot = mdot;   // stored for future steps (inertance or loss)
+                // ---- Determine air fraction for ghost ----
-
+                if (mdotFinal >= 0)
-                // Ghost state construction
+                {
-                float rhoFace = mdot >= 0 ? volRho : pipeRho;
-                float pFace = pFace0;
-                float uFace = MathF.Abs(mdot) / MathF.Max(rhoFace * area, 1e-12f);
-                float airFracGhost;
-                if (mdot >= 0)
                     airFracGhost = volAF;
+                }
                 else
                 {
                     airFracGhost = pipeAF;
-                    d.VolumePort.AirFraction = pipeAF;
+                    if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
                 }
 
-                if (mdot >= 0 && d.IsLeftEnd) uFace = +uFace;
+                // ---- Apply sign convention for velocity ----
-                else if (mdot >= 0 && !d.IsLeftEnd) uFace = -uFace;
+                if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace;
-                else if (mdot < 0 && d.IsLeftEnd) uFace = -uFace;
+                else if (mdotFinal >= 0 && !d.IsLeftEnd) uFace = -uFace;
-                else if (mdot < 0 && !d.IsLeftEnd) uFace = +uFace;
+                else if (mdotFinal < 0 && d.IsLeftEnd) uFace = -uFace;
+                else if (mdotFinal < 0 && !d.IsLeftEnd) uFace = +uFace;
 
+                // ---- Set ghost cells ----
                 if (d.IsLeftEnd)
                     _pipeSystem.SetGhostLeft(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
                 else
                     _pipeSystem.SetGhostRight(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
 
-                d.VolumePort.MassFlowRate = -mdot;
+                // ---- Update volume port (mass flow: positive into volume) ----
-                if (-mdot >= 0)
+                if (d.VolumePort != null)
                 {
-                    float pipeH = gamma / (gamma - 1f) * pipeP / MathF.Max(pipeRho, 1e-12f);
+                    d.VolumePort.MassFlowRate = -mdotFinal;
-                    d.VolumePort.SpecificEnthalpy = pipeH;
+
-                }
+                    // Set enthalpy of the stream entering the volume
-                else
+                    if (-mdotFinal >= 0)   // mass flowing into the volume (out of pipe)
-                {
+                    {
-                    d.VolumePort.SpecificEnthalpy = volH;
+                        float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
+                        d.VolumePort.SpecificEnthalpy = pipeH;
+                    }
+                    else   // mass flowing out of the volume (into pipe)
+                    {
+                        d.VolumePort.SpecificEnthalpy = volH;
+                    }
                 }
             }
         }
 
+        // ---------- Resolve open ends ----------
         public void ResolveOpenEnds(float dt)
         {
             for (int i = 0; i < OpenEndCount; i++)
@@ -303,6 +314,7 @@ namespace FluidSim.Core
 
                 bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
                 float rhoGhost, uGhost, pGhost, afGhost;
+
                 if (supersonic)
                 {
                     rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt;

Core/Pipesystem.cs

@@ -149,6 +149,7 @@ namespace FluidSim.Core
         public int GetPipeEnd(int pipeIdx) => _pipeEnd[pipeIdx];
         public float GetCellPressure(int i) => _p[i];
         public float GetCellDensity(int i) => _rho[i];
+        public float GetCellDx(int i) => _dx[i];
         public float GetCellVelocity(int i)
         {
             float rho = _rho[i];

Scenarios/HelmholtzScenario.cs

@@ -57,18 +57,23 @@ namespace FluidSim.Tests
             float rho0 = 101325f / (287f * 300f);
             pipeSystem = new PipeSystem(neckCells, pipeStart, pipeEnd, areas, dxs, rho0, 0f, 101325f);
 
-
-            // Energy loss
-            cavity.EnergyRelaxationRate = 80f;
-            pipeSystem.EnergyRelaxationRate = 0f;
-            pipeSystem.DampingMultiplier = 2000f;
-
             // --- Boundary system ---
             boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
 
+            float ComputeResistance(float decayTimeSeconds, float rho, float L_eff, float A)
+            {
+                // R = 2 * rho * L_eff / (A * decayTimeSeconds)
+                return 2f * rho * L_eff / (A * MathF.Max(decayTimeSeconds, 1e-6f));
+            }
+
             // Use steady orifice – the pipe already provides the inertia
-            boundaries.AddOrifice(cavityPort, pipeIndex: 0, isLeftEnd: true, areaIndex: cavityOrificeIdx, dischargeCoeff: 1f, lossCoefficient: 0.1f);
+            boundaries.AddOrificeWithInertance(
-            // LOSS COEFFICIENT BREAKS THE SYSTEM AT ~0.55, AT VALUES LOWER THAN THAT, IT SEEMS TO ONLY AFFECT VOLUME, NOT COMPOUND
+                cavityPort, pipeIndex: 0, isLeftEnd: true,
+                areaIndex: cavityOrificeIdx,
+                dischargeCoeff: 0.9f,
+                effectiveLength: neckLength,   // physical length (or L_eff)
+                lossCoefficient: 9000        // start with this, adjust for decay time
+            );
 
             // Open end at right side of pipe
             boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, neckArea);
@@ -78,7 +83,7 @@ namespace FluidSim.Tests
 
             // --- Solver ---
             // Slightly higher sub‑step count to ensure stability of the resonant oscillation
-            solver = new Solver { SubStepCount = 6, EnableProfiling = true };
+            solver = new Solver { SubStepCount = 6, EnableProfiling = false };
             solver.SetTimeStep(dt);
             solver.SetPipeSystem(pipeSystem);
             solver.SetBoundarySystem(boundaries);