FluidSim/Core/OrificeLink.cs

using System;
using FluidSim.Components;
using FluidSim.Interfaces;

namespace FluidSim.Core
{
    public class OrificeLink
    {
        public Port? VolumePort { get; }
        public Pipe1D Pipe { get; }
        public bool IsPipeLeftEnd { get; }
        public Func<double> AreaProvider { get; set; }
        public double DischargeCoefficient { get; set; } = 0.62;

        public double EffectiveLength { get; set; } = 0.001;
        public bool UseInertance { get; set; } = false;

        // Current mass flow (kg/s, positive = volume → pipe)
        private double _mdot;

        public double LastMassFlowRate { get; private set; }   // positive = into volume
        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)
        {
            VolumePort = volumePort;
            Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe));
            IsPipeLeftEnd = isPipeLeftEnd;
            AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider));
            _mdot = 0.0;
        }

        public void Resolve(double dtSub)
        {
            double area = AreaProvider();
            if (area < 1e-12 || VolumePort == null)
            {
                SetClosedWall();
                return;
            }

            // Gather states
            double volP = VolumePort.Pressure;
            double volRho = VolumePort.Density;
            double volT = VolumePort.Temperature;
            double volH = VolumePort.SpecificEnthalpy;

            (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd
                ? Pipe.GetInteriorStateLeft()
                : Pipe.GetInteriorStateRight();

            double pipeT = pipeP / Math.Max(pipeRho * 287.0, 1e-12);
            double gamma = 1.4;
            double R = 287.0;

            // ---- Steady‑state nozzle solution (gives correct exit state) ----
            double mdotSS;   // positive = volume → pipe
            double rhoFace0, uFace0, pFace0;
            if (volP >= pipeP)
            {
                IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, DischargeCoefficient,
                    out double mdotUpToDown, out rhoFace0, out uFace0, out pFace0);
                mdotSS = mdotUpToDown;
            }
            else
            {
                IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, DischargeCoefficient,
                    out double mdotUpToDown, out rhoFace0, out uFace0, out pFace0);
                mdotSS = -mdotUpToDown;
            }

            // ====== Hard physical cap: max sonic flow × 1.1 ======
            double upRho = mdotSS >= 0 ? volRho : pipeRho;
            double upT = mdotSS >= 0 ? volT : pipeT;
            double upC = Math.Sqrt(gamma * R * upT);
            double maxFlow = upRho * upC * area * 1.1;
            if (Math.Abs(mdotSS) > maxFlow)
                mdotSS = Math.Sign(mdotSS) * maxFlow;
            // ====================================================

            // ---- Dynamic update ----
            if (UseInertance)
            {
                double rhoUp = _mdot >= 0 ? volRho : pipeRho;
                double inertance = rhoUp * EffectiveLength / area;
                double dp = volP - pipeP;
                double resistance = Math.Abs(dp) / Math.Max(Math.Abs(mdotSS), 1e-12);
                double dmdot_dt = (dp - resistance * _mdot) / inertance;
                _mdot += dmdot_dt * dtSub;
            }
            else
            {
                _mdot = mdotSS;
            }

            // Clamp outflow to available mass (if finite volume)
            if (VolumePort.Owner is Volume0D vol)
            {
                double maxOut = vol.Mass / dtSub;
                if (_mdot > maxOut) _mdot = maxOut;
            }

            // ---- Ghost state ----
            double rhoFace = _mdot >= 0 ? volRho : pipeRho;
            double pFace = pFace0;
            double mdotMag = Math.Abs(_mdot);
            double uFace = mdotMag / (rhoFace * area);

            if (IsPipeLeftEnd)
                uFace = _mdot >= 0 ? uFace : -uFace;
            else
                uFace = _mdot >= 0 ? -uFace : uFace;

            if (IsPipeLeftEnd)
                Pipe.SetGhostLeft(rhoFace, uFace, pFace);
            else
                Pipe.SetGhostRight(rhoFace, uFace, pFace);

            // Store results (positive = into volume)
            LastMassFlowRate = -_mdot;
            LastFaceDensity = rhoFace;
            LastFaceVelocity = uFace;
            LastFacePressure = pFace;

            VolumePort.MassFlowRate = -_mdot;

            // Enthalpy transport
            if (-_mdot >= 0)   // inflow → pipe enthalpy
            {
                double hPipe = gamma / (gamma - 1.0) * pipeP / Math.Max(pipeRho, 1e-12);
                VolumePort.SpecificEnthalpy = hPipe;
            }
            else
            {
                VolumePort.SpecificEnthalpy = volH;
            }
        }

        private void SetClosedWall()
        {
            var (rInt, uInt, pInt) = IsPipeLeftEnd
                ? Pipe.GetInteriorStateLeft()
                : Pipe.GetInteriorStateRight();

            if (IsPipeLeftEnd)
                Pipe.SetGhostLeft(rInt, -uInt, pInt);
            else
                Pipe.SetGhostRight(rInt, -uInt, pInt);

            LastMassFlowRate = 0.0;
            LastFaceDensity = rInt;
            LastFaceVelocity = 0.0;
            LastFacePressure = pInt;
            if (VolumePort != null)
                VolumePort.MassFlowRate = 0.0;
        }
    }
}