FluidSim/Core/OrificeBoundary.cs

using System;
using FluidSim.Components;

namespace FluidSim.Core
{
    public static class OrificeBoundary
    {
        public static double MassFlow(double pA, double rhoA, double pB, double rhoB,
                                       Connection conn)
        {
            if (double.IsNaN(pA) || double.IsNaN(rhoA) || double.IsNaN(pB) || double.IsNaN(rhoB) ||
                double.IsInfinity(pA) || double.IsInfinity(rhoA) || double.IsInfinity(pB) || double.IsInfinity(rhoB) ||
                pA <= 0 || rhoA <= 0 || pB <= 0 || rhoB <= 0)
                return 0.0;

            double dp = pA - pB;
            double sign = Math.Sign(dp);
            double absDp = Math.Abs(dp);
            double rhoUp = dp >= 0 ? rhoA : rhoB;
            double pUp = dp >= 0 ? pA : pB;
            double pDown = dp >= 0 ? pB : pA;
            double delta = 1e-6 * pUp;

            if (absDp < delta)
            {
                double k = conn.DischargeCoefficient * conn.Area * Math.Sqrt(2 * rhoUp / delta);
                return k * dp;
            }
            else
            {
                double pr = pDown / pUp;
                double choked = Math.Pow(2.0 / (conn.Gamma + 1.0), conn.Gamma / (conn.Gamma - 1.0));
                if (pr < choked)
                {
                    double term = Math.Sqrt(conn.Gamma *
                        Math.Pow(2.0 / (conn.Gamma + 1.0), (conn.Gamma + 1.0) / (conn.Gamma - 1.0)));
                    double flow = conn.DischargeCoefficient * conn.Area *
                                  Math.Sqrt(rhoUp * pUp) * term;
                    return sign * flow;
                }
                else
                {
                    double ex = 1.0 - Math.Pow(pr, (conn.Gamma - 1.0) / conn.Gamma);
                    double flow = conn.DischargeCoefficient * conn.Area *
                        Math.Sqrt(2.0 * rhoUp * pUp * (conn.Gamma / (conn.Gamma - 1.0)) *
                                  pr * pr * ex);
                    return sign * flow;
                }
            }
        }

        public static void PipeVolumeFlux(double pPipe, double rhoPipe, double uPipe,
                                          double pVol, double rhoVol, double uVol,
                                          Connection conn, double pipeArea,
                                          bool isLeftBoundary,
                                          out double massFlux, out double momFlux, out double energyFlux)
        {
            // mass flow from pipe to volume (positive = pipe → volume)
            double mdot = MassFlow(pPipe, rhoPipe, pVol, rhoVol, conn);

            // Limit mass flow to the amount that can leave/enter the pipe cell
            double maxMdot = rhoPipe * pipeArea * 343.0;
            if (Math.Abs(mdot) > maxMdot) mdot = Math.Sign(mdot) * maxMdot;

            bool flowLeavesPipe = mdot > 0;
            double uFace, pFace, rhoFace;
            double massFluxPerArea;

            if (isLeftBoundary)
            {
                massFluxPerArea = -mdot / pipeArea;
                if (flowLeavesPipe)
                { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
                else
                { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
            }
            else // right boundary
            {
                massFluxPerArea = mdot / pipeArea;
                if (flowLeavesPipe)
                { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
                else
                { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
            }

            // Total enthalpy of the injected fluid (corrected: mass flux × total enthalpy)
            double specificEnthalpy = (1.4 / (1.4 - 1.0)) * pFace / Math.Max(rhoFace, 1e-12);
            double totalEnthalpy = specificEnthalpy + 0.5 * uFace * uFace;

            massFlux = massFluxPerArea;
            momFlux = massFluxPerArea * uFace + pFace;
            energyFlux = massFluxPerArea * totalEnthalpy;
        }
    }
}