using FluidSim.Components; using FluidSim.Interfaces; using System; 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; public int PipeIndex; public bool IsLeftEnd; public int AreaIndex; public float DischargeCoeff; // --- Inertance support --- public bool UseInertance; public float EffectiveLength; public float CurrentMdot; // kg/s, positive = volume → pipe // --- Loss coefficient (linear resistance) – inertance only --- // If 0 when UseInertance is true, a stable default is auto‑computed at runtime. public float LossCoefficient; // N·s/m⁵ or kg/(m⁴·s) } public struct OpenEndDesc { public int PipeIndex; public bool IsLeftEnd; public float AmbientPressure; public float Gamma; public float PipeArea; public float LastMassFlowRate; public float LastFacePressure; } private OrificeDesc[] _orifices; private OpenEndDesc[] _openEnds; private float[] _orificeAreas; private PipeSystem _pipeSystem; public BoundarySystem(PipeSystem pipeSystem, int maxOrifices, int maxOpenEnds) { _pipeSystem = pipeSystem; _orifices = new OrificeDesc[maxOrifices]; _openEnds = new OpenEndDesc[maxOpenEnds]; _orificeAreas = new float[maxOrifices]; } public int OrificeCount { get; private set; } public int OpenEndCount { get; private set; } // ---------- Add orifice (no inertance) ---------- // Simple isentropic nozzle – no built‑in loss. For dissipation use pipe damping // or the inertance model if you need a damped resonator. public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd, int areaIndex, float dischargeCoeff = 1f) { _orifices[OrificeCount] = new OrificeDesc { VolumePort = volumePort, PipeIndex = pipeIndex, IsLeftEnd = isLeftEnd, AreaIndex = areaIndex, DischargeCoeff = dischargeCoeff, UseInertance = false, EffectiveLength = 0f, CurrentMdot = 0f, LossCoefficient = 0f }; OrificeCount++; } // ---------- Add orifice with inertance ---------- // effectiveLength – length of the inertial slug (m), typically the physical neck length. // lossCoefficient – linear resistance (N·s/m⁵). If 0 (or omitted) an automatic stable // value will be computed from the pipe's characteristic impedance. public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd, int areaIndex, float dischargeCoeff, float effectiveLength, float lossCoefficient = 0f) { AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff); ref var d = ref _orifices[OrificeCount - 1]; d.UseInertance = true; d.EffectiveLength = effectiveLength; d.LossCoefficient = lossCoefficient; } public void AddOpenEnd(int pipeIndex, bool isLeftEnd, float ambientPressure, float pipeArea, float gamma = 1.4f) { int idx = OpenEndCount; _openEnds[idx] = new OpenEndDesc { PipeIndex = pipeIndex, IsLeftEnd = isLeftEnd, AmbientPressure = ambientPressure, Gamma = gamma, PipeArea = pipeArea }; OpenEndCount++; } public void SetOrificeAreas(float[] areas) { for (int i = 0; i < OrificeCount; i++) _orificeAreas[i] = areas[i]; } public float GetOpenEndMassFlow(int openEndIndex) { if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 0f; return _openEnds[openEndIndex].LastMassFlowRate; } public float GetOpenEndPressure(int openEndIndex) { if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 101325f; 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 as a wall ---- if (area < 1e-12f || d.VolumePort == null) { var (rInt, uInt, pInt) = d.IsLeftEnd ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex) : _pipeSystem.GetInteriorStateRight(d.PipeIndex); float afInt = d.IsLeftEnd ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex) : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex); if (d.IsLeftEnd) _pipeSystem.SetGhostLeft(d.PipeIndex, rInt, -uInt, pInt, afInt); else _pipeSystem.SetGhostRight(d.PipeIndex, rInt, -uInt, pInt, afInt); if (d.VolumePort != null) d.VolumePort.MassFlowRate = 0f; continue; } // ---- Preliminary isentropic solution ---- float mdotEst, rhoFaceEst, uFaceEst, pFaceEst; if (volP >= pipeP) { 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, Rgas, area, d.DischargeCoeff, out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst); mdotEst = -mdotEst; } // ---- Compute final mass flow with limiters ---- float mdotFinal, rhoFace, uFace, pFace, airFracGhost; if (d.UseInertance) { // ---- Inertance ODE with (possibly automatic) linear loss ---- float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho; float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f); float dp = volP - pipeP; float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f); float mdotNew = d.CurrentMdot + dmdot_dt * dt; // Limit outflow from volume (if volume owner is Volume0D) if (d.VolumePort.Owner is Volume0D vol0) { float maxOut = vol0.Mass / dt; if (mdotNew > maxOut) mdotNew = maxOut; if (mdotNew < -maxOut) mdotNew = -maxOut; } // Limit inflow from pipe – pipe cell cannot be emptied in one step { int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex) : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1; float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell); float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); // true cell area, not orifice float pipeDxAdj = _pipeSystem.GetCellDx(adjCell); float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj; float maxFromPipe = pipeCellMass / dt; if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe; } // Velocity clamp to Mach 0.9 float rhoFacePrelim = mdotNew >= 0 ? volRho : pipeRho; float uFacePrelim = MathF.Abs(mdotNew) / MathF.Max(rhoFacePrelim * area, 1e-12f); float cUp = mdotNew >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT); float maxU = 0.9f * cUp; if (uFacePrelim > maxU) { uFacePrelim = maxU; mdotNew = rhoFacePrelim * uFacePrelim * area * (mdotNew >= 0 ? 1f : -1f); } if (float.IsNaN(mdotNew)) mdotNew = 0f; d.CurrentMdot = mdotNew; mdotFinal = mdotNew; rhoFace = mdotFinal >= 0 ? volRho : pipeRho; pFace = pFaceEst; uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f); } else { // Standard quasi‑steady orifice mdotFinal = mdotEst; rhoFace = rhoFaceEst; uFace = uFaceEst; pFace = pFaceEst; // Limit outflow from volume (if Volume0D) if (d.VolumePort.Owner is Volume0D vol0) { float maxOut = vol0.Mass / dt; if (mdotFinal > maxOut) mdotFinal = maxOut; } // ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained ***** if (mdotFinal < 0) { int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex) : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1; float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell); float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); float pipeDxAdj = _pipeSystem.GetCellDx(adjCell); float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj; float maxFromPipe = pipeCellMass / dt; if (mdotFinal < -maxFromPipe) mdotFinal = -maxFromPipe; } d.CurrentMdot = mdotFinal; // Limit outflow from cylinder into pipe (positive mdot = volume → pipe) if (mdotFinal > 0f && d.VolumePort?.Owner is Cylinder cyl) { float maxOut = cyl.Mass / dt; if (mdotFinal > maxOut) mdotFinal = maxOut; } } // ---- Air fraction for ghost ---- if (mdotFinal >= 0) airFracGhost = volAF; else { airFracGhost = pipeAF; if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF; } // ---- 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); // ---- Update volume port ---- if (d.VolumePort != null) { d.VolumePort.MassFlowRate = -mdotFinal; if (-mdotFinal >= 0) // mass entering volume (out of pipe) { float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f); d.VolumePort.SpecificEnthalpy = pipeH; } else // mass leaving volume (into pipe) { d.VolumePort.SpecificEnthalpy = volH; } } } } // ---------- Resolve open ends ---------- public void ResolveOpenEnds(float dt) { for (int i = 0; i < OpenEndCount; i++) { ref var d = ref _openEnds[i]; var (rhoInt, uInt, pInt) = d.IsLeftEnd ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex) : _pipeSystem.GetInteriorStateRight(d.PipeIndex); float afInt = d.IsLeftEnd ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex) : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex); float gamma = d.Gamma; float gm1 = gamma - 1f; float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f)); float pAmb = d.AmbientPressure; // Characteristic solution (isentropic expansion to ambient) float Jplus = uInt + 2f * cInt / gm1; float Jminus = uInt - 2f * cInt / gm1; float s = pInt / MathF.Pow(rhoInt, gamma); float rhoIso = MathF.Pow(pAmb / s, 1f / gamma); float cIso = MathF.Sqrt(gamma * pAmb / MathF.Max(rhoIso, 1e-12f)); float uIso = d.IsLeftEnd ? (Jminus + 2f * cIso / gm1) : (Jplus - 2f * cIso / gm1); // Supersonic check bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt); if (!supersonic) { supersonic = d.IsLeftEnd ? (uIso <= -cIso) : (uIso >= cIso); } float rhoGhost, uGhost, pGhost, afGhost; if (supersonic) { rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt; } else { rhoGhost = rhoIso; uGhost = uIso; pGhost = pAmb; bool inflow = d.IsLeftEnd ? (uIso >= 0f) : (uIso <= 0f); afGhost = inflow ? 1f : afInt; } // ------- Mass flow limiter ------- int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex) : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1; float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell); float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); float pipeDxAdj = _pipeSystem.GetCellDx(adjCell); float cellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj; float area = d.PipeArea; float mdotRaw = rhoGhost * uGhost * area; // positive out of pipe if (d.IsLeftEnd) mdotRaw = -mdotRaw; // now positive = out of pipe // Outflow limit if (mdotRaw > 0 && mdotRaw * dt > cellMass) { mdotRaw = cellMass / dt; } // Inflow limit (allow up to 10× cell mass to avoid starving the pipe) else if (mdotRaw < 0 && -mdotRaw * dt > 10f * cellMass) { mdotRaw = -10f * cellMass / dt; } // Recompute uGhost from the limited mdot, keeping rhoGhost, pGhost float mdotMag = MathF.Abs(mdotRaw); uGhost = mdotMag / MathF.Max(rhoGhost * area, 1e-12f); if (d.IsLeftEnd) uGhost = (mdotRaw >= 0f) ? -uGhost : uGhost; else uGhost = (mdotRaw >= 0f) ? uGhost : -uGhost; // Apply ghost if (d.IsLeftEnd) _pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); else _pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); d.LastMassFlowRate = mdotRaw; d.LastFacePressure = pGhost; } } } }