From cf1bf30c8185e7c5442c31a078321935f82110b1 Mon Sep 17 00:00:00 2001 From: max Date: Sat, 9 May 2026 02:25:56 +0200 Subject: [PATCH] working helmholtz --- Core/BoundarySystem.cs | 200 +++++++++++++++++---------------- Core/Pipesystem.cs | 1 + Scenarios/HelmholtzScenario.cs | 23 ++-- 3 files changed, 121 insertions(+), 103 deletions(-) diff --git a/Core/BoundarySystem.cs b/Core/BoundarySystem.cs index 715d464..c66876c 100644 --- a/Core/BoundarySystem.cs +++ b/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 --- - public float LossCoefficient; // K factor for pressure drop = K * 0.5*rho*u^2 + // --- Loss coefficient (linear resistance) --- + 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 - 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 --- + // ---- Preliminary isentropic solution (used for face pressure if inertance is on) ---- 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 --- - float volP_eff = volP; - 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; + // ---- Compute ghost state ---- + float mdotFinal, rhoFace, uFace, pFace, airFracGhost; - // 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 dp = volP_eff - pipeP_eff; - float resistance = MathF.Abs(dp) / MathF.Max(MathF.Abs(mdotSS), 1e-12f); - float dmdot_dt = (dp - resistance * d.CurrentMdot) / inertance; - mdot = d.CurrentMdot + dmdot_dt * dt; + float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f); + float dp = volP - pipeP; + float Rlin = d.LossCoefficient; // linear resistance + // 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) - - // 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) + // ---- Determine air fraction for ghost ---- + if (mdotFinal >= 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; - else if (mdot >= 0 && !d.IsLeftEnd) uFace = -uFace; - else if (mdot < 0 && d.IsLeftEnd) uFace = -uFace; - else if (mdot < 0 && !d.IsLeftEnd) uFace = +uFace; + // ---- Apply sign convention for velocity ---- + if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace; + else if (mdotFinal >= 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; - if (-mdot >= 0) + // ---- Update volume port (mass flow: positive into volume) ---- + if (d.VolumePort != null) { - float pipeH = gamma / (gamma - 1f) * pipeP / MathF.Max(pipeRho, 1e-12f); - d.VolumePort.SpecificEnthalpy = pipeH; - } - else - { - d.VolumePort.SpecificEnthalpy = volH; + d.VolumePort.MassFlowRate = -mdotFinal; + + // Set enthalpy of the stream entering the volume + if (-mdotFinal >= 0) // mass flowing into the volume (out of pipe) + { + 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; diff --git a/Core/Pipesystem.cs b/Core/Pipesystem.cs index 349ebec..6a927e8 100644 --- a/Core/Pipesystem.cs +++ b/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]; diff --git a/Scenarios/HelmholtzScenario.cs b/Scenarios/HelmholtzScenario.cs index b62d2c7..18395c4 100644 --- a/Scenarios/HelmholtzScenario.cs +++ b/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); - // LOSS COEFFICIENT BREAKS THE SYSTEM AT ~0.55, AT VALUES LOWER THAN THAT, IT SEEMS TO ONLY AFFECT VOLUME, NOT COMPOUND + boundaries.AddOrificeWithInertance( + 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);