FluidSim/Components/Volume0D.cs

using System;

namespace FluidSim.Components
{
    public class Volume0D
    {
        public double Mass { get; set; }
        public double InternalEnergy { get; set; }

        public double Gamma { get; set; } = 1.4;
        public double GasConstant { get; set; } = 287.0;

        public double Volume { get; set; }
        public double Dvdt { get; set; }

        private double _dt;

        public double Density => Mass / Math.Max(Volume, 1e-12);
        public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12);
        public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12);

        public double MassFlowRateIn { get; set; }
        public double SpecificEnthalpyIn { get; set; }

        public Volume0D(double initialVolume, double initialPressure,
                        double initialTemperature, int sampleRate,
                        double gasConstant = 287.0, double gamma = 1.4)
        {
            GasConstant = gasConstant;
            Gamma = gamma;
            Volume = initialVolume;
            Dvdt = 0.0;
            _dt = 1.0 / sampleRate;

            double rho0 = initialPressure / (GasConstant * initialTemperature);
            Mass = rho0 * Volume;
            InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
        }

        public void Integrate(double dtOverride)
        {
            double dm = MassFlowRateIn * dtOverride;
            double dE = (MassFlowRateIn * SpecificEnthalpyIn) * dtOverride - Pressure * Dvdt * dtOverride;

            Mass += dm;
            InternalEnergy += dE;

            // ---- ABSOLUTE SAFEGUARD: keep at least 1 µg of gas at ambient pressure ----
            double minMass = 1e-9;
            double V = Math.Max(Volume, 1e-12);
            if (Mass < minMass)
            {
                Mass = minMass;
                InternalEnergy = 5000.0 * V / (Gamma - 1.0);   // 0.05 bar, not ambient
            }
            else if (InternalEnergy < 0.0)
            {
                InternalEnergy = 101325.0 * V / (Gamma - 1.0);
            }

            // Final non‑negative clamp
            if (Mass < 0.0) Mass = 0.0;
            if (InternalEnergy < 0.0) InternalEnergy = 0.0;
        }

        public void Integrate() => Integrate(_dt);
    }
}