FluidSim/Core/OpenEndLink.cs

using System;
using FluidSim.Components;

namespace FluidSim.Core
{
    /// <summary>
    /// Characteristic open‑end boundary condition after Jones (1978).
    /// For all subsonic flow (outflow and inflow), the ghost state is derived
    /// from the isentropic expansion to ambient pressure, using the pipe's entropy,
    /// and the outgoing Riemann invariant. This avoids a density jump at flow reversal.
    /// Supersonic outflow extrapolates the interior state.
    /// Now includes air fraction tracking: incoming air is fresh (AF=1), outgoing uses interior pipe AF.
    /// </summary>
    public class OpenEndLink
    {
        public Pipe1D Pipe { get; }
        public bool IsLeftEnd { get; }
        public double AmbientPressure { get; set; } = 101325.0;
        public double Gamma { get; set; } = 1.4;

        public double LastMassFlowRate { get; private set; }
        public double LastFaceDensity { get; private set; }
        public double LastFaceVelocity { get; private set; }
        public double LastFacePressure { get; private set; }

        public OpenEndLink(Pipe1D pipe, bool isLeftEnd)
        {
            Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe));
            IsLeftEnd = isLeftEnd;
        }

        public void Resolve(double dtSub)
        {
            (double rhoInt, double uInt, double pInt) = IsLeftEnd
                ? Pipe.GetInteriorStateLeft()
                : Pipe.GetInteriorStateRight();

            double airFracInt = IsLeftEnd
                ? Pipe.GetInteriorAirFractionLeft()
                : Pipe.GetInteriorAirFractionRight();

            double gamma = Gamma;
            double gm1 = gamma - 1.0;
            double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12));
            double pAmb = AmbientPressure;

            // Riemann invariants
            double J_plus = uInt + 2.0 * cInt / gm1;
            double J_minus = uInt - 2.0 * cInt / gm1;

            double rhoGhost, uGhost, pGhost, airFracGhost;

            // ---- Subsonic branch (used for both outflow and inflow) ----
            double s = pInt / Math.Pow(rhoInt, gamma);        // entropy constant
            double rhoIso = Math.Pow(pAmb / s, 1.0 / gamma);
            double cIso = Math.Sqrt(gamma * pAmb / Math.Max(rhoIso, 1e-12));
            double uIso = IsLeftEnd
                ? (J_minus + 2.0 * cIso / gm1)
                : (J_plus - 2.0 * cIso / gm1);

            // Check for supersonic outflow
            bool supersonic = IsLeftEnd
                ? (uInt <= -cInt)
                : (uInt >= cInt);

            if (!supersonic)
            {
                if (IsLeftEnd)
                    supersonic = uIso <= -cIso;
                else
                    supersonic = uIso >= cIso;
            }

            if (supersonic)
            {
                // Supersonic outflow – extrapolate interior
                rhoGhost = rhoInt;
                uGhost = uInt;
                pGhost = pInt;
                airFracGhost = airFracInt;   // whatever is leaving
            }
            else
            {
                // Subsonic flow – use isentropic state to ambient
                rhoGhost = rhoIso;
                uGhost = uIso;
                pGhost = pAmb;

                // Determine if inflow or outflow
                bool isInflow = IsLeftEnd ? (uIso >= 0) : (uIso <= 0); // positive u means into pipe from left end? Wait: left end: u>0 means flow to the right, into pipe. Right end: u>0 means flow to the right, out of pipe. Let's use mass flow sign later.
                // More straightforward: if using the isentropic state, the ghost velocity direction indicates flow. For inflow (ambient to pipe), airFraction = 1.0; for outflow, airFraction = interior's AF.
                if ((IsLeftEnd && uIso >= 0) || (!IsLeftEnd && uIso <= 0))
                {
                    // Inflow (ambient enters pipe)
                    airFracGhost = 1.0;
                }
                else
                {
                    // Outflow (pipe exits to ambient)
                    airFracGhost = airFracInt;
                }
            }

            // Apply ghost to pipe
            if (IsLeftEnd)
                Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost, airFracGhost);
            else
                Pipe.SetGhostRight(rhoGhost, uGhost, pGhost, airFracGhost);

            // Mass flow out of the pipe (positive = leaving)
            double area = Pipe.Area;
            double mdot = rhoGhost * uGhost * area;
            if (IsLeftEnd) mdot = -mdot;   // left end: positive u is into pipe, outward flow is -u
            LastMassFlowRate = mdot;
            LastFaceDensity = rhoGhost;
            LastFaceVelocity = uGhost;
            LastFacePressure = pGhost;
        }
    }
}