FluidSim/Components/Volume0D.cs

using System;
using System.Collections.Generic;
using FluidSim.Interfaces;

namespace FluidSim.Components
{
    public class Volume0D : IComponent
    {
        public List<Port> Ports { get; } = new List<Port>();

        private float _airMass;
        private float _exhaustMass;
        public float InternalEnergy;
        public float Volume;
        public float Dvdt;
        public float Gamma { get; set; } = 1.4f;
        public float GasConstant { get; set; } = 287f;
        public float AmbientPressure { get; set; } = 101325f;

        // ---------- Thermal relaxation to environment ----------
        /// <summary>Rate of heat transfer to the surroundings (1/s). 0 = adiabatic.</summary>
        public float EnergyRelaxationRate { get; set; } = 10f;
        /// <summary>Temperature to relax toward (K). Default is room temperature.</summary>
        public float AmbientTemperature { get; set; } = 300f;

        public float Mass => _airMass + _exhaustMass;
        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
        public float Density => Mass / MathF.Max(Volume, 1e-12f);
        public float Pressure => (Gamma - 1f) * InternalEnergy / MathF.Max(Volume, 1e-12f);
        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
        public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / MathF.Max(Density, 1e-12f);

        public Volume0D(float initialVolume, float initialPressure,
                        float initialTemperature, float gasConstant = 287f, float gamma = 1.4f)
        {
            GasConstant = gasConstant;
            Gamma = gamma;
            Volume = initialVolume;
            Dvdt = 0f;

            float rho0 = initialPressure / (GasConstant * initialTemperature);
            _airMass = rho0 * Volume;
            _exhaustMass = 0f;
            InternalEnergy = (initialPressure * Volume) / (Gamma - 1f);
        }

        public Port CreatePort()
        {
            var port = new Port { Owner = this };
            port.Pressure = Pressure;
            port.Density = Density;
            port.Temperature = Temperature;
            port.SpecificEnthalpy = SpecificEnthalpy;
            port.AirFraction = AirFraction;
            Ports.Add(port);
            return port;
        }

        public void SetPressure(float pressure, float? temperature = null)
        {
            float V = MathF.Max(Volume, 1e-12f);
            float T = temperature ?? Temperature;
            float rho = pressure / (GasConstant * T);
            float totalMass = rho * V;
            float af = AirFraction;
            _airMass = totalMass * af;
            _exhaustMass = totalMass * (1f - af);
            InternalEnergy = pressure * V / (Gamma - 1f);
        }

        public void UpdateState(float dt)
        {
            float totalMdotAir = 0f, totalMdotExhaust = 0f, totalEdot = 0f;
            foreach (var port in Ports)
            {
                float mdot = port.MassFlowRate;
                float af = mdot >= 0f ? port.AirFraction : AirFraction;
                totalMdotAir += mdot * af;
                totalMdotExhaust += mdot * (1f - af);
                totalEdot += mdot * port.SpecificEnthalpy;
            }

            float dAir = totalMdotAir * dt;
            float dExhaust = totalMdotExhaust * dt;
            float dE = totalEdot * dt - Pressure * Dvdt * dt;

            _airMass += dAir;
            _exhaustMass += dExhaust;
            InternalEnergy += dE;

            // ---- Thermal relaxation ----
            if (EnergyRelaxationRate > 0f)
            {
                float currentMass = Mass;
                if (currentMass > 1e-12f)
                {
                    // Target internal energy: current mass at ambient temperature
                    float targetE = currentMass * GasConstant * AmbientTemperature / (Gamma - 1f);
                    float relaxFactor = MathF.Exp(-EnergyRelaxationRate * dt);
                    InternalEnergy = targetE + (InternalEnergy - targetE) * relaxFactor;
                }
            }

            float V = MathF.Max(Volume, 1e-12f);
            float totalMass = _airMass + _exhaustMass;
            if (totalMass < 1e-9f)
            {
                _airMass = 1e-9f;
                _exhaustMass = 0f;
                InternalEnergy = AmbientPressure * V / (Gamma - 1f);
            }
            else if (InternalEnergy < 0f)
            {
                InternalEnergy = AmbientPressure * V / (Gamma - 1f);
            }

            if (_airMass < 0f) _airMass = 0f;
            if (_exhaustMass < 0f) _exhaustMass = 0f;

            float p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afr = AirFraction;
            foreach (var port in Ports)
            {
                port.Pressure = p;
                port.Density = rho;
                port.Temperature = T;
                port.SpecificEnthalpy = h;
                port.AirFraction = afr;
            }
        }

        IReadOnlyList<Port> IComponent.Ports => Ports;
    }
}