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 (for reference) ---- 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 ghost state ---- 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; // If loss coefficient not provided, use a tiny fraction of the pipe's characteristic impedance float Rlin = d.LossCoefficient; if (Rlin <= 0f) { // Auto‑sized linear drag: 0.5 % of Z_char float rhoRef = d.CurrentMdot >= 0 ? volRho : pipeRho; float cRef = d.CurrentMdot >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT); float Z_char = rhoRef * cRef / MathF.Max(area, 1e-12f); Rlin = 0.005f * Z_char; } float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f); float mdotNew = d.CurrentMdot + dmdot_dt * dt; // Symmetric flow limiters if (d.VolumePort.Owner is Volume0D vol0) { float maxOut = vol0.Mass / dt; if (mdotNew > maxOut) mdotNew = maxOut; if (mdotNew < -maxOut) mdotNew = -maxOut; } int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex) : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1; float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell); float pipeDxAdj = _pipeSystem.GetCellDx(adjCell); float pipeCellMass = pipeRhoAdj * area * pipeDxAdj; float maxFromPipe = pipeCellMass / dt; if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe; // Velocity clamp 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 (purely isentropic) ---- 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 (mdotFinal > maxOut) mdotFinal = maxOut; } // Safety velocity clamp (Mach 0.9) float cLocal = mdotFinal >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT); float maxULocal = 0.9f * cLocal; float uCheck = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f); if (uCheck > maxULocal) { uFace = maxULocal; mdotFinal = rhoFace * uFace * area * (mdotFinal >= 0 ? 1f : -1f); } d.CurrentMdot = mdotFinal; } // ---- Determine air fraction for ghost ---- if (mdotFinal >= 0) { airFracGhost = volAF; } else { airFracGhost = pipeAF; if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF; } // ---- 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); // ---- Update volume port ---- if (d.VolumePort != null) { d.VolumePort.MassFlowRate = -mdotFinal; if (-mdotFinal >= 0) // mass flowing into the volume { float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f); d.VolumePort.SpecificEnthalpy = pipeH; } else // mass flowing out of the volume { 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; 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); bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt); 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; } if (d.IsLeftEnd) _pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); else _pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); float area = d.PipeArea; float mdot = rhoGhost * uGhost * area; if (d.IsLeftEnd) mdot = -mdot; d.LastMassFlowRate = mdot; d.LastFacePressure = pGhost; } } } }