# Cryogenic LN2 Tank Module Implementation Plan > **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. **Goal:** Build a transient thermodynamic simulation of a cryogenic LN2 tank with helium pressurization, two-zone model (liquid + ullage), and plugin heat leak interface. **Architecture:** 3-variable ODE system [m_liq, U_liq, U_ull] solved by scipy.integrate.solve_ivp. Constant-pressure constraint provides He flow rate analytically. CoolProp for N2 properties, ideal gas for He. Independent module under src/cryo_tank/. **Tech Stack:** Python 3, CoolProp, scipy, numpy, matplotlib **Spec:** `docs/superpowers/specs/2026-04-16-cryo-tank-module-design.md` --- ## File Structure ``` src/cryo_tank/ __init__.py # Package marker (empty) config.py # All parameters: geometry, initial conditions, inlet/outlet, simulation control properties.py # CoolProp wrappers with AbstractState + lookup table + He ideal gas heat_leak.py # HeatLeakModel base + MLIHeatLeak + FoamHeatLeak tank_model.py # CryoTank class: geometry, state recovery, rhs(), He constraint solver.py # run() function: calls solve_ivp, returns history dict output.py # Plots (PNG) + history (NPZ) main.py # Entry point: assemble, run, output tests/cryo_tank/ __init__.py test_properties.py # CoolProp wrapper unit tests test_heat_leak.py # Heat leak model unit tests test_tank_model.py # Geometry, initial state, RHS unit tests test_integration.py # Short-run mass conservation + steady-state tests ``` --- ### Task 1: Package Scaffolding and Config **Files:** - Create: `src/cryo_tank/__init__.py` - Create: `src/cryo_tank/config.py` - Create: `tests/cryo_tank/__init__.py` - [ ] **Step 1: Create package directories and __init__ files** ```bash mkdir -p src/cryo_tank tests/cryo_tank touch src/cryo_tank/__init__.py tests/cryo_tank/__init__.py ``` - [ ] **Step 2: Write config.py** ```python # src/cryo_tank/config.py """ Configuration for the cryogenic LN2 tank simulation. Pure data module -- no functions, no side effects. """ import math # ---------- Gas constants ---------- R_UNIVERSAL = 8314.46 # J/(kmol*K) M_N2 = 28.014 # kg/kmol M_HE = 4.0026 # kg/kmol R_HE = R_UNIVERSAL / M_HE # 2077.1 J/(kg*K) R_N2 = R_UNIVERSAL / M_N2 # 296.8 J/(kg*K) # ---------- Tank geometry ---------- V_TOTAL = 420.1e-3 # m^3 (420.1 L) H_TANK = 0.5 # m (cylinder height) A_CROSS = V_TOTAL / H_TANK # m^2 (cross-section area) D_TANK = math.sqrt(4 * A_CROSS / math.pi) # m (diameter) A_SIDE = math.pi * D_TANK * H_TANK # m^2 (side wall) A_CAP = A_CROSS # m^2 (top or bottom cap) A_TOTAL = A_SIDE + 2 * A_CAP # m^2 (total surface) # ---------- Tank limits ---------- P_WORKING = 0.17e6 # Pa (working pressure, absolute) P_MAX = 0.8e6 # Pa (max bearing pressure) # ---------- Initial conditions ---------- T_INIT = 78.0 # K ULLAGE_FRACTION = 0.30 # gas pocket = 30% of V_TOTAL # ---------- Inlet / outlet ---------- MDOT_IN_LN2 = 1.144 # kg/s T_IN_LN2 = 77.0 # K MDOT_OUT_LN2 = 1.1895 # kg/s T_IN_HE = 100.0 # K # ---------- Heat transfer ---------- H_CONV_SURFACE = 50.0 # W/(m^2*K) liquid-to-ullage surface convection T_ENV = 300.0 # K ambient temperature # ---------- Simulation control ---------- T_END = 3600.0 # s (1 hour) RTOL = 1e-8 ATOL = 1e-10 # ---------- Output ---------- OUTPUT_DIR = "results/cryo_tank" # ---------- Validation ---------- assert V_TOTAL > 0 assert H_TANK > 0 assert 0 < ULLAGE_FRACTION < 1 assert P_WORKING > 0 assert P_MAX > P_WORKING assert T_INIT > 0 assert MDOT_IN_LN2 >= 0 assert MDOT_OUT_LN2 >= 0 assert T_IN_LN2 > 0 assert T_IN_HE > 0 assert H_CONV_SURFACE >= 0 assert T_ENV > 0 assert T_END > 0 ``` - [ ] **Step 3: Verify config imports cleanly** Run: `python3 -c "import sys; sys.path.insert(0,'src'); from cryo_tank.config import *; print(f'D={D_TANK:.3f}m, A_total={A_TOTAL:.3f}m2')"` Expected output: `D=1.034m, A_total=3.306m2` - [ ] **Step 4: Commit** ```bash git add src/cryo_tank/ tests/cryo_tank/ git commit -m "feat(cryo_tank): package scaffolding and config module" ``` --- ### Task 2: CoolProp Property Wrappers (properties.py) **Files:** - Create: `src/cryo_tank/properties.py` - Create: `tests/cryo_tank/test_properties.py` - [ ] **Step 1: Write failing tests for properties** ```python # tests/cryo_tank/test_properties.py """Tests for CoolProp property wrappers.""" import pytest import sys sys.path.insert(0, "src") from cryo_tank.properties import ( ln2_rho, ln2_h, ln2_u, ln2_T_from_u, n2_vapor_u, n2_sat_pressure, he_u, he_h, he_cp, he_cv, ) from cryo_tank.config import P_WORKING class TestLN2Properties: """Liquid nitrogen properties at P = 0.17 MPa.""" def test_ln2_density_at_78K(self): rho = ln2_rho(78.0, P_WORKING) assert 800 < rho < 810 # ~803 kg/m3 def test_ln2_enthalpy_at_77K(self): h = ln2_h(77.0, P_WORKING) assert -130000 < h < -110000 # ~-122695 J/kg def test_ln2_internal_energy_at_78K(self): u = ln2_u(78.0, P_WORKING) assert -130000 < u < -110000 # ~-120865 J/kg def test_ln2_T_from_u_roundtrip(self): T_orig = 78.0 u = ln2_u(T_orig, P_WORKING) T_recovered = ln2_T_from_u(u, P_WORKING) assert abs(T_recovered - T_orig) < 0.01 class TestN2VaporProperties: """N2 vapor properties.""" def test_n2_sat_pressure_at_78K(self): P_sat = n2_sat_pressure(78.0) assert 0.10e6 < P_sat < 0.12e6 # ~0.1093 MPa def test_n2_vapor_internal_energy_at_78K(self): u = n2_vapor_u(78.0) assert 50000 < u < 60000 # ~55547 J/kg class TestHeliumProperties: """Helium (ideal gas) properties.""" def test_he_cp_near_5196(self): cp = he_cp() assert abs(cp - 5196.2) < 10 # monatomic ideal gas def test_he_cv_near_3117(self): cv = he_cv() assert abs(cv - 3117.1) < 10 def test_he_enthalpy_at_100K(self): h = he_h(100.0) # He enthalpy ~ cp * T (relative to some reference) # CoolProp gives ~524762 J/kg at 100K assert 500000 < h < 550000 def test_he_internal_energy_at_78K(self): u = he_u(78.0) assert 200000 < u < 280000 # ~247932 J/kg ``` - [ ] **Step 2: Run tests to verify they fail** Run: `python3 -m pytest tests/cryo_tank/test_properties.py -v` Expected: FAIL (module not found) - [ ] **Step 3: Implement properties.py** ```python # src/cryo_tank/properties.py """ Fluid property wrappers for liquid nitrogen, N2 vapor, and helium. Performance strategy (per spec Section 8.1): - N2 (liquid & vapor): CoolProp with persistent AbstractState objects - He: analytical ideal gas (cp=5196.2 J/(kg*K), cv=3117.1 J/(kg*K)) - Lookup tables for the ODE hot path (built at import time) """ import numpy as np import CoolProp.CoolProp as CP from CoolProp import AbstractState # --------------------------------------------------------------------------- # Persistent CoolProp AbstractState objects (reused across calls) # --------------------------------------------------------------------------- _n2_state = AbstractState("HEOS", "Nitrogen") _he_state = AbstractState("HEOS", "Helium") # --------------------------------------------------------------------------- # Liquid nitrogen (LN2) properties at a given (T, P) # --------------------------------------------------------------------------- def ln2_rho(T, P): """LN2 density [kg/m^3].""" _n2_state.update(CP.PT_INPUTS, P, T) return _n2_state.rhomass() def ln2_h(T, P): """LN2 specific enthalpy [J/kg].""" _n2_state.update(CP.PT_INPUTS, P, T) return _n2_state.hmass() def ln2_u(T, P): """LN2 specific internal energy [J/kg].""" _n2_state.update(CP.PT_INPUTS, P, T) return _n2_state.umass() def ln2_T_from_u(u, P): """Recover LN2 temperature from specific internal energy [K]. Uses lookup table interpolation for speed; falls back to CoolProp if outside the table range. """ return float(np.interp(u, _ln2_u_table, _ln2_T_table)) # --------------------------------------------------------------------------- # N2 vapor properties (at saturation or specified conditions) # --------------------------------------------------------------------------- def n2_sat_pressure(T): """N2 saturation pressure [Pa] at temperature T.""" _n2_state.update(CP.QT_INPUTS, 1.0, T) return _n2_state.p() def n2_vapor_u(T): """N2 saturated vapor specific internal energy [J/kg] at temperature T.""" _n2_state.update(CP.QT_INPUTS, 1.0, T) return _n2_state.umass() def n2_vapor_rho(T): """N2 saturated vapor density [kg/m^3] at temperature T.""" _n2_state.update(CP.QT_INPUTS, 1.0, T) return _n2_state.rhomass() # --------------------------------------------------------------------------- # Helium properties (ideal gas: cp=5/2 R, cv=3/2 R, monatomic) # --------------------------------------------------------------------------- _HE_CP = 5196.2 # J/(kg*K), = 5/2 * R_He _HE_CV = 3117.1 # J/(kg*K), = 3/2 * R_He # Reference state: CoolProp He at T_ref=0K gives u_ref, h_ref # We match CoolProp's reference by computing offset at a known point. _he_state.update(CP.PT_INPUTS, 170000.0, 100.0) _HE_H_REF = _he_state.hmass() - _HE_CP * 100.0 # h = cp*T + h_ref _HE_U_REF = _he_state.umass() - _HE_CV * 100.0 # u = cv*T + u_ref def he_cp(): """He specific heat at constant pressure [J/(kg*K)].""" return _HE_CP def he_cv(): """He specific heat at constant volume [J/(kg*K)].""" return _HE_CV def he_h(T): """He specific enthalpy [J/kg] (ideal gas).""" return _HE_CP * T + _HE_H_REF def he_u(T): """He specific internal energy [J/kg] (ideal gas).""" return _HE_CV * T + _HE_U_REF def he_T_from_u(u): """Recover He temperature from specific internal energy [K].""" return (u - _HE_U_REF) / _HE_CV # --------------------------------------------------------------------------- # Lookup table for LN2: u(T) -> T at P = 0.17 MPa (built at import time) # --------------------------------------------------------------------------- _LN2_T_MIN = 65.0 _LN2_T_MAX = 82.0 # stay below saturation at 0.17 MPa (~82.03 K) _LN2_TABLE_N = 200 _P_WORK = 170000.0 _ln2_T_table = np.linspace(_LN2_T_MIN, _LN2_T_MAX, _LN2_TABLE_N) _ln2_u_table = np.array([ln2_u(T, _P_WORK) for T in _ln2_T_table]) # _ln2_u_table is monotonically increasing, so np.interp works for inverse lookup ``` - [ ] **Step 4: Run tests to verify they pass** Run: `python3 -m pytest tests/cryo_tank/test_properties.py -v` Expected: All 8 tests PASS - [ ] **Step 5: Commit** ```bash git add src/cryo_tank/properties.py tests/cryo_tank/test_properties.py git commit -m "feat(cryo_tank): CoolProp property wrappers with He ideal gas and LN2 lookup table" ``` --- ### Task 3: Heat Leak Models (heat_leak.py) **Files:** - Create: `src/cryo_tank/heat_leak.py` - Create: `tests/cryo_tank/test_heat_leak.py` - [ ] **Step 1: Write failing tests** ```python # tests/cryo_tank/test_heat_leak.py """Tests for heat leak models.""" import sys sys.path.insert(0, "src") from cryo_tank.heat_leak import HeatLeakModel, MLIHeatLeak, FoamHeatLeak class TestMLIHeatLeak: def test_mli_returns_constant_heat_flux(self): model = MLIHeatLeak(A_total=3.306, q_mli=1.5) Q = model.compute(T_inner=78.0, T_env=300.0) assert abs(Q - 3.306 * 1.5) < 1e-10 def test_mli_default_q_is_1(self): model = MLIHeatLeak(A_total=3.306) Q = model.compute(T_inner=78.0, T_env=300.0) assert abs(Q - 3.306) < 1e-10 def test_mli_independent_of_temperature(self): model = MLIHeatLeak(A_total=3.306, q_mli=2.0) Q1 = model.compute(T_inner=78.0, T_env=300.0) Q2 = model.compute(T_inner=80.0, T_env=250.0) assert abs(Q1 - Q2) < 1e-10 class TestFoamHeatLeak: def test_foam_constant_k(self): model = FoamHeatLeak(A_total=3.306, k_eff=0.03, delta=0.05) Q = model.compute(T_inner=78.0, T_env=300.0) expected = 3.306 * 0.03 * (300.0 - 78.0) / 0.05 assert abs(Q - expected) < 1e-6 def test_foam_callable_k(self): def k_func(T): return 0.01 + 0.0001 * T # linear k(T) model = FoamHeatLeak(A_total=3.306, k_eff=k_func, delta=0.05) Q = model.compute(T_inner=78.0, T_env=300.0) T_mean = (300.0 + 78.0) / 2.0 k_at_mean = k_func(T_mean) expected = 3.306 * k_at_mean * (300.0 - 78.0) / 0.05 assert abs(Q - expected) < 1e-6 def test_foam_zero_dT_gives_zero_Q(self): model = FoamHeatLeak(A_total=3.306, k_eff=0.03, delta=0.05) Q = model.compute(T_inner=300.0, T_env=300.0) assert abs(Q) < 1e-10 ``` - [ ] **Step 2: Run tests to verify they fail** Run: `python3 -m pytest tests/cryo_tank/test_heat_leak.py -v` Expected: FAIL - [ ] **Step 3: Implement heat_leak.py** ```python # src/cryo_tank/heat_leak.py """ Heat leak models for the cryogenic tank. Provides a plugin interface (HeatLeakModel base class) and two built-in implementations: MLI (vacuum multi-layer) and Foam (wrap insulation). """ class HeatLeakModel: """Base class for heat leak models. Subclasses must implement compute(T_inner, T_env) -> Q [W]. Positive Q means heat flows INTO the tank. """ def compute(self, T_inner, T_env): raise NotImplementedError class MLIHeatLeak(HeatLeakModel): """Vacuum multi-layer insulation. Heat flux is approximately constant (independent of temperature) in the typical cryogenic operating range. Parameters ---------- A_total : float Total tank surface area [m^2]. q_mli : float, default 1.0 Specific heat flux [W/m^2]. """ def __init__(self, A_total, q_mli=1.0): self.A_total = A_total self.q_mli = q_mli def compute(self, T_inner, T_env): return self.A_total * self.q_mli class FoamHeatLeak(HeatLeakModel): """Foam or wrap insulation with 1D steady conduction model. Parameters ---------- A_total : float Total tank surface area [m^2]. k_eff : float or callable Effective thermal conductivity [W/(m*K)]. If callable, signature k_eff(T) -> float, evaluated at T_mean. delta : float Insulation thickness [m]. """ def __init__(self, A_total, k_eff, delta): self.A_total = A_total self._k_eff = k_eff self.delta = delta def _get_k(self, T_mean): if callable(self._k_eff): return self._k_eff(T_mean) return self._k_eff def compute(self, T_inner, T_env): T_mean = (T_inner + T_env) / 2.0 k = self._get_k(T_mean) return self.A_total * k * (T_env - T_inner) / self.delta ``` - [ ] **Step 4: Run tests to verify they pass** Run: `python3 -m pytest tests/cryo_tank/test_heat_leak.py -v` Expected: All 6 tests PASS - [ ] **Step 5: Commit** ```bash git add src/cryo_tank/heat_leak.py tests/cryo_tank/test_heat_leak.py git commit -m "feat(cryo_tank): heat leak models (MLI + Foam with k(T) support)" ``` --- ### Task 4: Tank Model -- Geometry and Initialization (tank_model.py) **Files:** - Create: `src/cryo_tank/tank_model.py` - Create: `tests/cryo_tank/test_tank_model.py` - [ ] **Step 1: Write failing tests for geometry and initial state** ```python # tests/cryo_tank/test_tank_model.py """Tests for CryoTank model.""" import pytest import sys sys.path.insert(0, "src") import numpy as np from cryo_tank.tank_model import CryoTank from cryo_tank.heat_leak import MLIHeatLeak from cryo_tank.config import ( V_TOTAL, H_TANK, A_CROSS, A_TOTAL, P_WORKING, T_INIT, ULLAGE_FRACTION, MDOT_IN_LN2, T_IN_LN2, MDOT_OUT_LN2, T_IN_HE, H_CONV_SURFACE, T_ENV, ) def _make_tank(): """Create a CryoTank with default config and MLI heat leak.""" heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0) return CryoTank( V_total=V_TOTAL, H_tank=H_TANK, P_work=P_WORKING, T_init=T_INIT, ullage_fraction=ULLAGE_FRACTION, mdot_in_ln2=MDOT_IN_LN2, T_in_ln2=T_IN_LN2, mdot_out_ln2=MDOT_OUT_LN2, T_in_he=T_IN_HE, h_conv=H_CONV_SURFACE, T_env=T_ENV, heat_leak_model=heat_leak, ) class TestGeometry: def test_cross_section_area(self): tank = _make_tank() assert abs(tank.A_cross - 0.8402) < 0.001 def test_total_surface_area(self): tank = _make_tank() assert abs(tank.A_total - 3.306) < 0.01 def test_wetted_area_at_70_percent_fill(self): tank = _make_tank() level = 0.7 * H_TANK # 0.35 m A_wet, A_dry = tank.wetted_areas(level) # A_wet = bottom cap + side * level expected_wet = A_CROSS + np.pi * tank.D * level assert abs(A_wet - expected_wet) < 0.01 assert abs(A_wet + A_dry - A_TOTAL) < 0.01 class TestInitialState: def test_initial_liquid_mass(self): tank = _make_tank() y0 = tank.initial_state() m_liq = y0[0] # rho_LN2(78K, 0.17MPa) ~ 803.3 kg/m3, V_liq = 0.2941 m3 assert 235 < m_liq < 237 # ~236.25 kg def test_initial_fill_fraction(self): tank = _make_tank() y0 = tank.initial_state() info = tank.derive(y0) assert abs(info['fill_fraction'] - 0.70) < 0.01 def test_initial_temperatures(self): tank = _make_tank() y0 = tank.initial_state() info = tank.derive(y0) assert abs(info['T_liq'] - T_INIT) < 0.1 assert abs(info['T_ull'] - T_INIT) < 1.0 def test_initial_pressure_components_sum_to_P_working(self): tank = _make_tank() y0 = tank.initial_state() info = tank.derive(y0) P_N2 = info['P_N2'] P_He = info['P_He'] assert abs(P_N2 + P_He - P_WORKING) / P_WORKING < 1e-6 ``` - [ ] **Step 2: Run tests to verify they fail** Run: `python3 -m pytest tests/cryo_tank/test_tank_model.py -v` Expected: FAIL - [ ] **Step 3: Implement tank_model.py (geometry + initialization + derive)** ```python # src/cryo_tank/tank_model.py """ CryoTank: two-zone (liquid + ullage) cryogenic tank model. State vector y = [m_liq, U_liq, U_ull] (3 components). Derived quantities (T, V, m_He, etc.) computed by derive(y). ODE right-hand side provided by rhs(t, y). """ import math import warnings import numpy as np from cryo_tank.config import R_HE, R_N2 from cryo_tank import properties as prop class CryoTank: """Two-zone cryogenic LN2 tank with He pressurization. Parameters ---------- V_total : float Total tank volume [m^3] H_tank : float Cylinder height [m] P_work : float Working pressure [Pa] T_init : float Initial temperature [K] (both zones) ullage_fraction : float Initial gas volume / total volume mdot_in_ln2 : float LN2 inlet mass flow [kg/s] T_in_ln2 : float LN2 inlet temperature [K] mdot_out_ln2 : float LN2 outlet mass flow [kg/s] T_in_he : float He inlet temperature [K] h_conv : float Surface heat transfer coeff [W/(m^2*K)] T_env : float Environment temperature [K] heat_leak_model : HeatLeakModel Plugin for heat leak calculation """ def __init__(self, V_total, H_tank, P_work, T_init, ullage_fraction, mdot_in_ln2, T_in_ln2, mdot_out_ln2, T_in_he, h_conv, T_env, heat_leak_model): # Geometry self.V_total = V_total self.H_tank = H_tank self.A_cross = V_total / H_tank self.D = math.sqrt(4 * self.A_cross / math.pi) self.A_side = math.pi * self.D * H_tank self.A_cap = self.A_cross self.A_total = self.A_side + 2 * self.A_cap # Operating conditions self.P_work = P_work self.mdot_in_ln2 = mdot_in_ln2 self.T_in_ln2 = T_in_ln2 self.mdot_out_ln2 = mdot_out_ln2 self.T_in_he = T_in_he self.h_conv = h_conv self.T_env = T_env self.heat_leak_model = heat_leak_model # Precompute constant inlet enthalpies self.h_in_ln2 = prop.ln2_h(T_in_ln2, P_work) self.h_in_he = prop.he_h(T_in_he) # Net liquid flow (constant) self.dm_liq_dt = mdot_in_ln2 - mdot_out_ln2 # Initial state computation self._T_init = T_init self._ullage_fraction = ullage_fraction V_liq_0 = (1.0 - ullage_fraction) * V_total V_ull_0 = ullage_fraction * V_total # Liquid initial state rho_liq_0 = prop.ln2_rho(T_init, P_work) self._m_liq_0 = rho_liq_0 * V_liq_0 self._U_liq_0 = self._m_liq_0 * prop.ln2_u(T_init, P_work) # Ullage initial state: N2 vapor at saturation + He to fill pressure P_N2_0 = prop.n2_sat_pressure(T_init) self.m_N2_ull = prop.n2_vapor_rho(T_init) * V_ull_0 # FIXED for all time P_He_0 = P_work - P_N2_0 self._m_He_0 = P_He_0 * V_ull_0 / (R_HE * T_init) U_N2_ull_0 = self.m_N2_ull * prop.n2_vapor_u(T_init) U_He_0 = self._m_He_0 * prop.he_u(T_init) self._U_ull_0 = U_N2_ull_0 + U_He_0 # Saturation temperature warning threshold self._T_sat = 82.0 # approximate, from CoolProp: ~82.03 K at 0.17 MPa def initial_state(self): """Return the ODE initial state vector y0 = [m_liq, U_liq, U_ull].""" return np.array([self._m_liq_0, self._U_liq_0, self._U_ull_0]) def wetted_areas(self, liquid_level): """Return (A_wet, A_dry) for the given liquid level [m].""" level = max(0.0, min(liquid_level, self.H_tank)) A_wet = self.A_cap + math.pi * self.D * level A_dry = self.A_cap + math.pi * self.D * (self.H_tank - level) return A_wet, A_dry def derive(self, y): """Compute all derived quantities from state vector y. Returns a dict with T_liq, T_ull, V_liq, V_ull, liquid_level, fill_fraction, m_He, P_N2, P_He, etc. """ m_liq, U_liq, U_ull = y[0], y[1], y[2] # Liquid zone u_liq = U_liq / m_liq # specific internal energy T_liq = prop.ln2_T_from_u(u_liq, self.P_work) rho_liq = prop.ln2_rho(T_liq, self.P_work) V_liq = m_liq / rho_liq liquid_level = V_liq / self.A_cross fill_fraction = liquid_level / self.H_tank # Ullage zone V_ull = self.V_total - V_liq # N2 partial pressure (ideal gas for vapor in ullage) P_N2 = self.m_N2_ull * R_N2 * 78.0 / V_ull # initial approx # Better: iterate to find T_ull first, then P_N2 # For now, solve T_ull from ullage energy # He mass from pressure constraint # First estimate T_ull from U_ull assuming P_N2 ~ const P_He = self.P_work - P_N2 # Ullage internal energy: U_ull = m_N2_ull * u_N2(T_ull) + m_He * u_He(T_ull) # For N2 vapor, u_N2 ~ cv_N2 * T_ull (ideal gas approx) # For He, u_He = cv_He * T_ull + u_ref # This is implicit in T_ull. Use iterative approach: # Start with T_ull guess, compute m_He, recompute T_ull from energy. T_ull = self._solve_ullage_temperature(U_ull, V_ull) # Recompute P_N2 and m_He with correct T_ull P_N2 = self.m_N2_ull * R_N2 * T_ull / V_ull P_He = self.P_work - P_N2 m_He = P_He * V_ull / (R_HE * T_ull) return { 'T_liq': T_liq, 'T_ull': T_ull, 'V_liq': V_liq, 'V_ull': V_ull, 'liquid_level': liquid_level, 'fill_fraction': fill_fraction, 'rho_liq': rho_liq, 'm_He': m_He, 'P_N2': P_N2, 'P_He': P_He, } def _solve_ullage_temperature(self, U_ull, V_ull): """Solve for T_ull given total ullage internal energy and volume. U_ull = m_N2_ull * u_N2_vap(T) + m_He(T) * u_He(T) where m_He(T) = (P_work - m_N2_ull * R_N2 * T / V_ull) * V_ull / (R_He * T) Solved by Newton iteration. """ T = self._T_init # initial guess for _ in range(50): P_N2 = self.m_N2_ull * R_N2 * T / V_ull P_He = self.P_work - P_N2 if P_He < 0: P_He = 0.0 m_He = P_He * V_ull / (R_HE * T) # N2 vapor internal energy (ideal gas approx): u = cv_N2 * T + u_ref # Use CoolProp reference: u_N2_vap(78K) = 55546.6, cv_N2 ~ 743 J/(kg*K) cv_N2 = 743.0 u_N2_ref = prop.n2_vapor_u(78.0) - cv_N2 * 78.0 u_N2 = cv_N2 * T + u_N2_ref U_calc = self.m_N2_ull * u_N2 + m_He * prop.he_u(T) residual = U_calc - U_ull if abs(residual) < 1e-3: # converged (< 1 mJ) return T # Numerical derivative dT = 0.01 P_N2_p = self.m_N2_ull * R_N2 * (T + dT) / V_ull P_He_p = max(0.0, self.P_work - P_N2_p) m_He_p = P_He_p * V_ull / (R_HE * (T + dT)) u_N2_p = cv_N2 * (T + dT) + u_N2_ref U_calc_p = self.m_N2_ull * u_N2_p + m_He_p * prop.he_u(T + dT) dU_dT = (U_calc_p - U_calc) / dT if abs(dU_dT) < 1e-20: break T = T - residual / dU_dT T = max(50.0, min(T, 400.0)) # clamp warnings.warn(f"_solve_ullage_temperature did not converge: T={T:.2f}, residual={residual:.2e}") return T def rhs(self, t, y): """ODE right-hand side: dy/dt = [dm_liq/dt, dU_liq/dt, dU_ull/dt]. This is called by scipy.integrate.solve_ivp. """ info = self.derive(y) m_liq = y[0] T_liq = info['T_liq'] T_ull = info['T_ull'] V_ull = info['V_ull'] rho_liq = info['rho_liq'] liquid_level = info['liquid_level'] m_He = info['m_He'] # --- Heat transfer --- Q_liq_to_ull = self.h_conv * self.A_cross * (T_liq - T_ull) Q_leak = self.heat_leak_model.compute(T_liq, self.T_env) A_wet, A_dry = self.wetted_areas(liquid_level) A_total = A_wet + A_dry Q_leak_liq = Q_leak * A_wet / A_total if A_total > 0 else 0.0 Q_leak_ull = Q_leak * A_dry / A_total if A_total > 0 else 0.0 # --- He flow rate (analytical, per spec Section 2.6) --- mdot_He = self._solve_he_flow_rate(info, Q_liq_to_ull, Q_leak_ull) # --- Liquid zone --- dm_liq_dt = self.dm_liq_dt # constant: mdot_in - mdot_out h_liq = prop.ln2_h(T_liq, self.P_work) dU_liq_dt = (self.mdot_in_ln2 * self.h_in_ln2 - self.mdot_out_ln2 * h_liq - Q_liq_to_ull + Q_leak_liq) # --- Ullage zone --- dU_ull_dt = mdot_He * self.h_in_he + Q_liq_to_ull + Q_leak_ull # --- Warnings --- if T_liq > self._T_sat - 1.0: warnings.warn( f"T_liq={T_liq:.2f}K approaching saturation ({self._T_sat:.1f}K); " "evaporation effects may be significant." ) return np.array([dm_liq_dt, dU_liq_dt, dU_ull_dt]) def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull): """Analytically solve for m_dot_He from dP/dt = 0 constraint. The key equation: P_total = P_N2 + P_He = const. P_N2 = m_N2_ull * R_N2 * T_ull / V_ull (N2 ideal gas, m_N2_ull = const) P_He = m_He * R_HE * T_ull / V_ull dP/dt = 0 implies dP_He/dt = -dP_N2/dt. Expanding and solving for m_dot_He gives a linear equation. See spec Section 2.6 for full derivation. """ T_ull = info['T_ull'] V_ull = info['V_ull'] m_He = info['m_He'] rho_liq = info['rho_liq'] # dV_ull/dt = -dV_liq/dt = -dm_liq/dt / rho_liq = -(mdot_in - mdot_out) / rho_liq dV_ull_dt = -self.dm_liq_dt / rho_liq # positive when liquid drains # Total ullage heat input (excluding He inlet, which we're solving for) Q_ull_no_he = Q_liq_to_ull + Q_leak_ull # Ullage total cv*mass (for dT_ull/dt estimation) cv_N2 = 743.0 # J/(kg*K), N2 vapor cv_He = prop.he_cv() C_ull = self.m_N2_ull * cv_N2 + m_He * cv_He # total heat capacity [J/K] # From dP_total/dt = 0 and ideal gas for both species: # P_total * dV_ull/dt + V_ull * dP_total/dt = d(P_total * V_ull)/dt # Since dP_total/dt = 0: # We need: d(n_total * R_u * T_ull)/dt = P_total * dV_ull/dt # where n_total = m_N2/M_N2 + m_He/M_He (in moles) # # Simplified linear solve: # m_dot_He = [P_work * dV_ull/dt - (m_N2*R_N2 + m_He*R_HE) * dT_ull_no_he / C_ull * V_ull/T_ull] # / [R_HE * T_ull - h_in_he * (m_N2*R_N2 + m_He*R_HE) * V_ull / (C_ull * T_ull)] # # More directly: from P_total = (m_N2*R_N2 + m_He*R_HE) * T_ull / V_ull # dP/dt = 0 = R_HE*T_ull/V_ull * dm_He/dt # + (m_N2*R_N2 + m_He*R_HE)/V_ull * dT_ull/dt # - (m_N2*R_N2 + m_He*R_HE)*T_ull/V_ull^2 * dV_ull/dt # # dT_ull/dt = (m_dot_He * h_in_he + Q_ull_no_he) / C_ull (from energy eq) # (note: this is an approximation; strictly dT = dU/C with PdV work, but # for the constraint solve it's sufficient) R_mix = self.m_N2_ull * R_N2 + m_He * R_HE # effective "mR" [J/K] # Coefficient of m_dot_He in dP/dt = 0: # A * m_dot_He + B = 0 # m_dot_He = -B / A A = R_HE * T_ull / V_ull + R_mix / (V_ull * C_ull) * self.h_in_he B = R_mix / (V_ull * C_ull) * Q_ull_no_he - R_mix * T_ull / (V_ull ** 2) * dV_ull_dt if abs(A) < 1e-30: return 0.0 mdot_He = -B / A # Clamp: He can only flow in (strict mode per spec Section 10.3) if mdot_He < 0: warnings.warn( f"He backflow requested (m_dot_He={mdot_He:.4e} kg/s); " "clamping to 0. Pressure may drift above target." ) mdot_He = 0.0 return mdot_He ``` - [ ] **Step 4: Run tests to verify they pass** Run: `python3 -m pytest tests/cryo_tank/test_tank_model.py -v` Expected: All 5 tests PASS - [ ] **Step 5: Commit** ```bash git add src/cryo_tank/tank_model.py tests/cryo_tank/test_tank_model.py git commit -m "feat(cryo_tank): CryoTank model with geometry, initialization, derive, rhs, He constraint" ``` --- ### Task 5: ODE Solver Driver (solver.py) **Files:** - Create: `src/cryo_tank/solver.py` - [ ] **Step 1: Implement solver.py** ```python # src/cryo_tank/solver.py """ ODE solver driver for the cryogenic tank simulation. Calls scipy.integrate.solve_ivp with the CryoTank.rhs method. Returns a history dict with all output quantities as time series. """ import warnings import numpy as np from scipy.integrate import solve_ivp def _liquid_empty_event(t, y): """Event function: triggers when m_liq reaches 0.""" return y[0] # m_liq _liquid_empty_event.terminal = True _liquid_empty_event.direction = -1 def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0): """Run the cryogenic tank simulation. Parameters ---------- tank : CryoTank Configured tank model instance. t_end : float End time [s]. rtol, atol : float ODE solver tolerances. max_step : float Maximum time step [s]. Returns ------- dict History with keys: 't', 'T_liq', 'T_ull', 'm_liq', 'm_He', 'mdot_He', 'V_liq', 'V_ull', 'liquid_level', 'fill_fraction', 'Q_leak', 'Q_leak_liq', 'Q_leak_ull', 'Q_liq_to_ull'. """ y0 = tank.initial_state() sol = solve_ivp( tank.rhs, [0.0, t_end], y0, method='RK45', rtol=rtol, atol=atol, max_step=max_step, events=[_liquid_empty_event], dense_output=True, ) if not sol.success: raise RuntimeError(f"ODE solver failed: {sol.message}") if sol.t_events[0].size > 0: warnings.warn( f"Tank emptied at t = {sol.t_events[0][0]:.1f} s " f"(before t_end = {t_end:.1f} s)" ) # --- Post-process: compute derived quantities at each output time --- t = sol.t n = len(t) history = { 't': t, 'm_liq': sol.y[0], 'U_liq': sol.y[1], 'U_ull': sol.y[2], 'T_liq': np.zeros(n), 'T_ull': np.zeros(n), 'm_He': np.zeros(n), 'mdot_He': np.zeros(n), 'V_liq': np.zeros(n), 'V_ull': np.zeros(n), 'liquid_level': np.zeros(n), 'fill_fraction': np.zeros(n), 'Q_leak': np.zeros(n), 'Q_leak_liq': np.zeros(n), 'Q_leak_ull': np.zeros(n), 'Q_liq_to_ull': np.zeros(n), } for i in range(n): y_i = sol.y[:, i] info = tank.derive(y_i) history['T_liq'][i] = info['T_liq'] history['T_ull'][i] = info['T_ull'] history['m_He'][i] = info['m_He'] history['V_liq'][i] = info['V_liq'] history['V_ull'][i] = info['V_ull'] history['liquid_level'][i] = info['liquid_level'] history['fill_fraction'][i] = info['fill_fraction'] # Recompute heat terms for recording T_liq = info['T_liq'] T_ull = info['T_ull'] level = info['liquid_level'] Q_liq_to_ull = tank.h_conv * tank.A_cross * (T_liq - T_ull) Q_leak = tank.heat_leak_model.compute(T_liq, tank.T_env) A_wet, A_dry = tank.wetted_areas(level) A_sum = A_wet + A_dry Q_leak_liq = Q_leak * A_wet / A_sum if A_sum > 0 else 0.0 Q_leak_ull = Q_leak * A_dry / A_sum if A_sum > 0 else 0.0 history['Q_liq_to_ull'][i] = Q_liq_to_ull history['Q_leak'][i] = Q_leak history['Q_leak_liq'][i] = Q_leak_liq history['Q_leak_ull'][i] = Q_leak_ull # He flow rate via finite difference on m_He (post-processing only, not in RHS) dt = np.diff(t) dm_He = np.diff(history['m_He']) history['mdot_He'][0] = dm_He[0] / dt[0] if len(dt) > 0 else 0.0 history['mdot_He'][1:] = dm_He / dt return history ``` - [ ] **Step 2: Commit** ```bash git add src/cryo_tank/solver.py git commit -m "feat(cryo_tank): ODE solver driver with liquid-empty event detection" ``` --- ### Task 6: Integration Tests **Files:** - Create: `tests/cryo_tank/test_integration.py` - [ ] **Step 1: Write integration tests** ```python # tests/cryo_tank/test_integration.py """Integration tests for the cryogenic tank simulation.""" import sys sys.path.insert(0, "src") import numpy as np from cryo_tank.tank_model import CryoTank from cryo_tank.heat_leak import MLIHeatLeak from cryo_tank.solver import run from cryo_tank.config import ( V_TOTAL, H_TANK, P_WORKING, T_INIT, ULLAGE_FRACTION, MDOT_IN_LN2, T_IN_LN2, MDOT_OUT_LN2, T_IN_HE, H_CONV_SURFACE, T_ENV, A_TOTAL, ) def _make_tank(**overrides): """Create a tank with default config, allowing overrides.""" kw = dict( V_total=V_TOTAL, H_tank=H_TANK, P_work=P_WORKING, T_init=T_INIT, ullage_fraction=ULLAGE_FRACTION, mdot_in_ln2=MDOT_IN_LN2, T_in_ln2=T_IN_LN2, mdot_out_ln2=MDOT_OUT_LN2, T_in_he=T_IN_HE, h_conv=H_CONV_SURFACE, T_env=T_ENV, heat_leak_model=MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0), ) kw.update(overrides) return CryoTank(**kw) class TestMassConservation: def test_liquid_mass_change_matches_net_flow(self): """Over a short run, dm_liq should equal (mdot_in - mdot_out) * dt.""" tank = _make_tank() history = run(tank, t_end=10.0, max_step=1.0) m_liq_0 = history['m_liq'][0] m_liq_f = history['m_liq'][-1] t_f = history['t'][-1] expected_dm = (MDOT_IN_LN2 - MDOT_OUT_LN2) * t_f actual_dm = m_liq_f - m_liq_0 rel_err = abs(actual_dm - expected_dm) / abs(expected_dm) assert rel_err < 1e-6, f"Mass conservation error: rel_err={rel_err:.2e}" class TestSteadyState: def test_zero_flow_zero_leak_is_static(self): """With no flow and no heat leak, state should not change.""" tank = _make_tank( mdot_in_ln2=0.0, mdot_out_ln2=0.0, h_conv=0.0, heat_leak_model=MLIHeatLeak(A_total=A_TOTAL, q_mli=0.0), ) history = run(tank, t_end=100.0, max_step=10.0) T_liq = history['T_liq'] T_ull = history['T_ull'] assert abs(T_liq[-1] - T_liq[0]) < 0.01, f"T_liq drifted: {T_liq[0]:.3f} -> {T_liq[-1]:.3f}" assert abs(T_ull[-1] - T_ull[0]) < 0.1, f"T_ull drifted: {T_ull[0]:.3f} -> {T_ull[-1]:.3f}" class TestPhysicalBehavior: def test_liquid_level_decreases(self): """With net outflow, liquid level should decrease.""" tank = _make_tank() history = run(tank, t_end=60.0, max_step=5.0) assert history['fill_fraction'][-1] < history['fill_fraction'][0] def test_he_flow_rate_positive(self): """He should always flow in (pressurization), not out.""" tank = _make_tank() history = run(tank, t_end=60.0, max_step=5.0) assert np.all(history['mdot_He'] >= -1e-10) # allow tiny numerical noise ``` - [ ] **Step 2: Run all tests** Run: `python3 -m pytest tests/cryo_tank/ -v` Expected: All tests PASS (properties: 8, heat_leak: 6, tank_model: 5, integration: 4 = 23 total) - [ ] **Step 3: Commit** ```bash git add tests/cryo_tank/test_integration.py git commit -m "test(cryo_tank): integration tests for mass conservation, steady state, physical behavior" ``` --- ### Task 7: Output and Visualization (output.py) **Files:** - Create: `src/cryo_tank/output.py` - [ ] **Step 1: Implement output.py** ```python # src/cryo_tank/output.py """ Output helpers for the cryogenic tank simulation. Generates PNG plots and NPZ data files. """ import os import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt def save_history(history, path): """Save all history time series to a compressed .npz file.""" dirname = os.path.dirname(path) if dirname: os.makedirs(dirname, exist_ok=True) np.savez_compressed(path, **history) def plot_temperatures(history, path): """Plot T_liq and T_ull vs time.""" fig, ax = plt.subplots(figsize=(10, 5)) t = history['t'] ax.plot(t, history['T_liq'], label='T_liq (liquid)') ax.plot(t, history['T_ull'], label='T_ull (ullage)') ax.set_xlabel('Time [s]') ax.set_ylabel('Temperature [K]') ax.set_title('Tank temperatures vs time') ax.grid(True) ax.legend() fig.tight_layout() fig.savefig(path, dpi=120) plt.close(fig) def plot_liquid_level(history, path): """Plot fill fraction and liquid level vs time.""" fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True) t = history['t'] ax1.plot(t, history['fill_fraction'] * 100) ax1.set_ylabel('Fill fraction [%]') ax1.set_title('Liquid level vs time') ax1.grid(True) ax2.plot(t, history['liquid_level'] * 1000) ax2.set_xlabel('Time [s]') ax2.set_ylabel('Liquid level [mm]') ax2.grid(True) fig.tight_layout() fig.savefig(path, dpi=120) plt.close(fig) def plot_he_flow(history, path): """Plot helium mass and flow rate vs time.""" fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True) t = history['t'] ax1.plot(t, history['m_He'] * 1000) ax1.set_ylabel('He mass [g]') ax1.set_title('Helium pressurization vs time') ax1.grid(True) ax2.plot(t, history['mdot_He'] * 1000) ax2.set_xlabel('Time [s]') ax2.set_ylabel('He flow rate [g/s]') ax2.grid(True) fig.tight_layout() fig.savefig(path, dpi=120) plt.close(fig) def plot_heat_fluxes(history, path): """Plot all heat transfer terms vs time.""" fig, ax = plt.subplots(figsize=(10, 5)) t = history['t'] ax.plot(t, history['Q_leak'], label='Q_leak (total)') ax.plot(t, history['Q_leak_liq'], label='Q_leak_liq', linestyle='--') ax.plot(t, history['Q_leak_ull'], label='Q_leak_ull', linestyle='--') ax.plot(t, history['Q_liq_to_ull'], label='Q_liq_to_ull') ax.set_xlabel('Time [s]') ax.set_ylabel('Heat flux [W]') ax.set_title('Heat transfer vs time') ax.grid(True) ax.legend() fig.tight_layout() fig.savefig(path, dpi=120) plt.close(fig) ``` - [ ] **Step 2: Commit** ```bash git add src/cryo_tank/output.py git commit -m "feat(cryo_tank): output module with temperature, level, He flow, and heat flux plots" ``` --- ### Task 8: Entry Point (main.py) and End-to-End Run **Files:** - Create: `src/cryo_tank/main.py` - [ ] **Step 1: Implement main.py** ```python # src/cryo_tank/main.py """ Entry point for the cryogenic LN2 tank simulation. Run from project root: python3 src/cryo_tank/main.py """ import os import sys _HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.dirname(_HERE)) from cryo_tank.config import ( V_TOTAL, H_TANK, P_WORKING, T_INIT, ULLAGE_FRACTION, MDOT_IN_LN2, T_IN_LN2, MDOT_OUT_LN2, T_IN_HE, H_CONV_SURFACE, T_ENV, A_TOTAL, T_END, RTOL, ATOL, OUTPUT_DIR, ) from cryo_tank.tank_model import CryoTank from cryo_tank.heat_leak import MLIHeatLeak from cryo_tank.solver import run from cryo_tank.output import ( save_history, plot_temperatures, plot_liquid_level, plot_he_flow, plot_heat_fluxes, ) def main(): os.makedirs(OUTPUT_DIR, exist_ok=True) # --- Assemble --- heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0) tank = CryoTank( V_total=V_TOTAL, H_tank=H_TANK, P_work=P_WORKING, T_init=T_INIT, ullage_fraction=ULLAGE_FRACTION, mdot_in_ln2=MDOT_IN_LN2, T_in_ln2=T_IN_LN2, mdot_out_ln2=MDOT_OUT_LN2, T_in_he=T_IN_HE, h_conv=H_CONV_SURFACE, T_env=T_ENV, heat_leak_model=heat_leak, ) y0 = tank.initial_state() info0 = tank.derive(y0) print(f"Initial state:") print(f" m_liq = {y0[0]:.2f} kg") print(f" T_liq = {info0['T_liq']:.2f} K, T_ull = {info0['T_ull']:.2f} K") print(f" fill_fraction = {info0['fill_fraction']:.1%}") print(f" m_He = {info0['m_He']*1000:.2f} g") print(f" P_N2 = {info0['P_N2']/1e6:.4f} MPa, P_He = {info0['P_He']/1e6:.4f} MPa") print(f"Running to t_end = {T_END:.0f} s ...") print() # --- Run --- history = run(tank, t_end=T_END, rtol=RTOL, atol=ATOL) n_steps = len(history['t']) t_final = history['t'][-1] print(f"Simulation complete: {n_steps} output points, t_final = {t_final:.1f} s") print(f" T_liq: {history['T_liq'][0]:.2f} -> {history['T_liq'][-1]:.2f} K") print(f" T_ull: {history['T_ull'][0]:.2f} -> {history['T_ull'][-1]:.2f} K") print(f" fill_fraction: {history['fill_fraction'][0]:.1%} -> {history['fill_fraction'][-1]:.1%}") print(f" m_He: {history['m_He'][0]*1000:.2f} -> {history['m_He'][-1]*1000:.2f} g") print(f" T_out (LN2 outlet) = T_liq = {history['T_liq'][-1]:.2f} K") # --- Output --- save_history(history, os.path.join(OUTPUT_DIR, "cryo_tank_history.npz")) plot_temperatures(history, os.path.join(OUTPUT_DIR, "cryo_tank_temperatures.png")) plot_liquid_level(history, os.path.join(OUTPUT_DIR, "cryo_tank_level.png")) plot_he_flow(history, os.path.join(OUTPUT_DIR, "cryo_tank_he_flow.png")) plot_heat_fluxes(history, os.path.join(OUTPUT_DIR, "cryo_tank_heat.png")) print(f"\nOutputs written to {OUTPUT_DIR}/") if __name__ == "__main__": main() ``` - [ ] **Step 2: Run the full simulation** Run: `python3 src/cryo_tank/main.py` Expected output (approximate): ``` Initial state: m_liq = 236.25 kg T_liq = 78.00 K, T_ull = 78.00 K fill_fraction = 70.0% m_He = 47.24 g P_N2 = 0.1093 MPa, P_He = 0.0607 MPa Running to t_end = 3600 s ... Simulation complete: ... output points, t_final = 3600.0 s T_liq: 78.00 -> ... K T_ull: 78.00 -> ... K fill_fraction: 70.0% -> ...% ... ``` - [ ] **Step 3: Verify output files exist** Run: `ls -lh results/cryo_tank/` Expected: 5 files (1 npz + 4 png) - [ ] **Step 4: Run full test suite** Run: `python3 -m pytest tests/ -v` Expected: All tests pass (pipe system tests + cryo_tank tests) - [ ] **Step 5: Commit** ```bash git add src/cryo_tank/main.py git commit -m "feat(cryo_tank): entry point with full simulation pipeline" ```