diff --git a/src/gas_cylinder.py b/src/gas_cylinder.py new file mode 100644 index 0000000..d73bea4 --- /dev/null +++ b/src/gas_cylinder.py @@ -0,0 +1,171 @@ +""" +High-pressure gas cylinder model with configurable working fluid. + +The cylinder is a 0D lumped model. Its primary state is total mass and +total internal energy; pressure, temperature, density, enthalpy, and other +thermodynamic quantities are recovered from CoolProp on demand. +""" + +import numpy as np +import CoolProp.CoolProp as CP +from CoolProp import AbstractState + + +DEFAULT_FLUID = "Helium" +DEFAULT_P_INIT = 18.031e6 # Pa, absolute +DEFAULT_VOLUME = 20e-3 # m^3, 20 L +DEFAULT_T_INIT = 78.0 # K + + +class HighPressureGasCylinder: + """0D high-pressure gas cylinder. + + Parameters + ---------- + fluid : str + CoolProp fluid name, for example ``"Helium"`` or ``"Nitrogen"``. + P_init : float + Initial absolute pressure [Pa]. + V : float + Cylinder volume [m^3]. + T_init : float + Initial temperature [K]. + backend : str + CoolProp backend name. ``"HEOS"`` is used by default. + """ + + def __init__(self, fluid=DEFAULT_FLUID, P_init=DEFAULT_P_INIT, + V=DEFAULT_VOLUME, T_init=DEFAULT_T_INIT, backend="HEOS"): + if not fluid: + raise ValueError("fluid must be a non-empty CoolProp fluid name") + if P_init <= 0: + raise ValueError("P_init must be > 0") + if V <= 0: + raise ValueError("V must be > 0") + if T_init <= 0: + raise ValueError("T_init must be > 0") + + self.fluid = fluid + self.backend = backend + self.V = V + self._state = AbstractState(backend, fluid) + + self._state.update(CP.PT_INPUTS, P_init, T_init) + rho = self._state.rhomass() + self.mass = rho * V + self.U = self.mass * self._state.umass() + + def _update_state(self): + rho = self.rho + u = self.specific_internal_energy + self._state.update(CP.DmassUmass_INPUTS, rho, u) + return self._state + + @property + def rho(self): + """Gas density [kg/m^3].""" + return self.mass / self.V + + @property + def specific_internal_energy(self): + """Specific internal energy [J/kg].""" + return self.U / self.mass + + @property + def P(self): + """Absolute pressure [Pa].""" + return self._update_state().p() + + @property + def T(self): + """Temperature [K].""" + return self._update_state().T() + + @property + def h(self): + """Specific enthalpy [J/kg].""" + return self._update_state().hmass() + + @property + def cp(self): + """Specific heat at constant pressure [J/(kg K)].""" + return self._update_state().cpmass() + + @property + def cv(self): + """Specific heat at constant volume [J/(kg K)].""" + return self._update_state().cvmass() + + @property + def gamma(self): + """Local heat capacity ratio cp/cv [-].""" + return self.cp / self.cv + + @property + def R_specific(self): + """Specific gas constant [J/(kg K)].""" + state = self._update_state() + return state.gas_constant() / state.molar_mass() + + @property + def compressibility_factor(self): + """Compressibility factor Z = P/(rho R T) [-].""" + return self.P / (self.rho * self.R_specific * self.T) + + def ghost_state(self): + """Return a stagnant conservative state ``[rho, rho*u, rho*E]``. + + This is useful as a boundary ghost cell for conservative flow models. + The state uses the real-fluid internal energy density. A pipe solver + that assumes an ideal-gas EOS must still be checked for EOS consistency + before coupling it directly to this real-fluid cylinder. + """ + return np.array([self.rho, 0.0, self.rho * self.specific_internal_energy]) + + def apply_flux(self, mdot, edot, dt, sign): + """Update cylinder mass and energy from a boundary flux. + + Parameters + ---------- + mdot : float + Mass flow rate [kg/s]. + edot : float + Energy flow rate [W]. + dt : float + Time step [s]. + sign : int + ``-1`` for outflow from the cylinder, ``+1`` for inflow. + """ + if dt < 0: + raise ValueError("dt must be >= 0") + if sign not in (-1, 1): + raise ValueError("sign must be +1 or -1") + + self.mass += sign * mdot * dt + self.U += sign * edot * dt + + if self.mass <= 0: + raise RuntimeError( + f"Cylinder mass non-positive after apply_flux: mass={self.mass}, " + f"mdot={mdot}, edot={edot}, dt={dt}, sign={sign}" + ) + + # Force a thermodynamic validity check immediately after the update. + self._update_state() + + def state_summary(self): + """Return common cylinder state quantities as a dictionary.""" + return { + "fluid": self.fluid, + "V": self.V, + "mass": self.mass, + "rho": self.rho, + "P": self.P, + "T": self.T, + "U": self.U, + "u": self.specific_internal_energy, + "h": self.h, + "gamma": self.gamma, + "R_specific": self.R_specific, + "Z": self.compressibility_factor, + } diff --git a/tests/test_gas_cylinder.py b/tests/test_gas_cylinder.py new file mode 100644 index 0000000..469097a --- /dev/null +++ b/tests/test_gas_cylinder.py @@ -0,0 +1,75 @@ +import numpy as np +import pytest + +from gas_cylinder import ( + DEFAULT_FLUID, + DEFAULT_P_INIT, + DEFAULT_T_INIT, + DEFAULT_VOLUME, + HighPressureGasCylinder, +) + + +def test_default_high_pressure_cylinder_state_matches_inputs(): + cylinder = HighPressureGasCylinder() + + assert cylinder.fluid == DEFAULT_FLUID + assert cylinder.V == DEFAULT_VOLUME + assert abs(cylinder.P - DEFAULT_P_INIT) / DEFAULT_P_INIT < 1e-8 + assert abs(cylinder.T - DEFAULT_T_INIT) < 1e-8 + assert cylinder.mass > 0 + assert cylinder.U > 0 + assert cylinder.gamma > 1.0 + assert cylinder.compressibility_factor > 0.0 + + +def test_cylinder_accepts_configurable_fluid_pressure_volume_temperature(): + cylinder = HighPressureGasCylinder( + fluid="Nitrogen", + P_init=2.0e6, + V=50e-3, + T_init=300.0, + ) + + assert cylinder.fluid == "Nitrogen" + assert abs(cylinder.P - 2.0e6) / 2.0e6 < 1e-8 + assert abs(cylinder.T - 300.0) < 1e-8 + assert abs(cylinder.V - 50e-3) < 1e-15 + assert cylinder.mass > 0 + + +def test_cylinder_ghost_state_is_stagnant_conservative_vector(): + cylinder = HighPressureGasCylinder() + ghost = cylinder.ghost_state() + + assert ghost.shape == (3,) + assert abs(ghost[0] - cylinder.rho) < 1e-12 + assert ghost[1] == 0.0 + assert abs(ghost[2] - cylinder.rho * cylinder.specific_internal_energy) < 1e-6 + + +def test_cylinder_apply_flux_outflow_reduces_mass_energy_and_pressure(): + cylinder = HighPressureGasCylinder() + mass_before = cylinder.mass + U_before = cylinder.U + P_before = cylinder.P + + mdot = 1e-3 + edot = mdot * cylinder.h + dt = 1.0 + cylinder.apply_flux(mdot=mdot, edot=edot, dt=dt, sign=-1) + + assert abs(cylinder.mass - (mass_before - mdot * dt)) < 1e-12 + assert abs(cylinder.U - (U_before - edot * dt)) < 1e-6 + assert cylinder.P < P_before + + +def test_cylinder_rejects_invalid_inputs(): + with pytest.raises(ValueError): + HighPressureGasCylinder(fluid="") + with pytest.raises(ValueError): + HighPressureGasCylinder(P_init=0.0) + with pytest.raises(ValueError): + HighPressureGasCylinder(V=0.0) + with pytest.raises(ValueError): + HighPressureGasCylinder(T_init=0.0)