lujingze/pipe-system-simulation-test
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
master
pipe-system-simulation-test/docs/superpowers/plans/2026-04-16-cryo-tank-module.md
T
ljz 0f57ba94f3 代码仓库移植
2026-06-03 15:41:04 +08:00

1433 lines
44 KiB
Markdown
Raw Permalink Blame History

# 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"
```
Reference in new issue View Git Blame Copy Permalink