pipe-system-simulation-test/docs/superpowers/specs/2026-04-16-cryo-tank-module-design.md

Cryogenic LN2 Tank Simulation Module -- Design Specification

Date: 2026-04-16 Status: Approved

1. Overview

A transient thermodynamic simulation module for a cryogenic liquid nitrogen (LN2) storage tank with helium pressurization. The module is an independent Python package within the existing pipe-system-simulation project, designed for future coupling with the pipe system.

1.1 Physical Scenario

A cylindrical cryogenic tank stores liquid nitrogen. It has:

• Two inlets: LN2 inlet and He pressurization inlet
• One outlet: LN2 outlet
• A heat leak interface for connecting insulation models

The tank pressure is maintained at a constant 0.17 MPa (absolute) by adjusting helium flow. Liquid nitrogen flows in at 1.144 kg/s (77 K) and out at 1.1895 kg/s. The helium inlet temperature is 100 K, and its flow rate is determined by the constant-pressure constraint.

1.2 Key Design Decisions (from brainstorming)

Decision	Choice	Rationale
Relationship to pipe system	Independent module (A)	Can run standalone, future coupling possible
Tank model	Two-zone (liquid + ullage) (A)	Captures liquid/gas temperature difference
Interphase mass transfer	Not considered	Only heat transfer between liquid and gas zones
Fluid properties	CoolProp (C)	Highest accuracy for cryogenic conditions
Heat leak interface	Plugin-style (C)	Flexible, extensible, decoupled
Pressure control	Algebraic constraint (A)	Exact constant pressure, no PID tuning
Time integration	scipy solve_ivp (A)	Adaptive stepping, good accuracy
Liquid temperature	Uniform (homogeneous)	Justified by continuous flow and shallow liquid layer

2. Physical Model

2.1 Two-Zone Model

The tank is divided into:

• Liquid zone (lower): LN2 at uniform temperature T_liq
• Ullage zone (upper): mixture of N2 vapor + He gas at uniform temperature T_ull

The two zones exchange heat across the liquid surface. No evaporation or condensation (no interphase mass transfer).

2.2 State Variables (ODE)

The ODE state vector has 3 components: y = [m_liq, U_liq, U_ull]

• m_liq [kg]: liquid nitrogen mass
• U_liq [J]: liquid zone total internal energy
• U_ull [J]: ullage zone total internal energy (N2 vapor + He combined)

2.3 Derived Quantities (not ODE states)

These are computed at each evaluation from the state + constraints:

• T_liq: liquid temperature (from m_liq, U_liq via CoolProp)
• T_ull: ullage temperature (from U_ull, m_N2_ull, m_He via CoolProp/ideal gas)
• V_liq = m_liq / rho_LN2(T_liq, P): liquid volume
• V_ull = V_total - V_liq: ullage volume
• liquid_level = V_liq / A_cross: liquid height in the cylinder
• m_He: helium mass from constant-pressure constraint (see Section 2.6)
• m_dot_He: helium mass flow rate (time derivative of m_He, also from constraint)
• T_out = T_liq: LN2 outlet temperature (uniform assumption)

2.4 Liquid Zone Governing Equations

Mass conservation:

dm_liq/dt = m_dot_in_LN2 - m_dot_out_LN2
          = 1.144 - 1.1895
          = -0.0455 kg/s  (constant)

Energy conservation:

dU_liq/dt = m_dot_in_LN2 * h_in_LN2
          - m_dot_out_LN2 * h_liq
          - Q_liq_to_ull
          + Q_leak_liq

Where:

• h_in_LN2 = h_N2_liquid(T_in=77K, P=0.17MPa) from CoolProp [J/kg]
• h_liq = h_N2_liquid(T_liq, P=0.17MPa) from CoolProp [J/kg] -- also the outlet enthalpy
• Q_liq_to_ull: heat transfer from liquid to ullage across the liquid surface [W]
• Q_leak_liq: heat leak into liquid zone from environment [W]

2.5 Ullage Zone Governing Equations

Mass conservation:

m_N2_ull = constant (no mass transfer, set at initialization)
dm_He/dt = m_dot_He (determined by constant-pressure constraint)

Energy conservation:

dU_ull/dt = m_dot_He * h_in_He
          + Q_liq_to_ull
          + Q_leak_ull

Where:

• h_in_He = h_He(T_in=100K, P=0.17MPa) from CoolProp [J/kg]
• Q_liq_to_ull: heat transfer from liquid surface (same magnitude, opposite sign as in liquid equation)
• Q_leak_ull: heat leak into ullage zone from environment [W]

2.6 Constant-Pressure Constraint (Analytical m_dot_He Derivation)

Tank pressure is maintained at P_total = 0.17 MPa at all times. The ullage gas follows Dalton's law:

P_total = P_N2 + P_He = 0.17 MPa

The N2 partial pressure depends on the fixed N2 vapor mass, ullage volume, and ullage temperature. The He mass required to provide the remaining pressure:

P_N2 = f(m_N2_ull, V_ull, T_ull)  (CoolProp or ideal gas)
P_He = P_total - P_N2
m_He = P_He * V_ull / (R_He * T_ull)  (He is well-approximated as ideal gas at these conditions)

Where R_He = R_universal / M_He = 8314.46 / 4.0026 = 2077.1 J/(kg*K).

m_He is NOT an ODE state variable. It is a derived quantity from the algebraic constraint.

Analytical derivation of m_dot_He (avoiding finite-difference instability):

Since P_total = const, differentiating dP/dt = 0 and using the ideal gas relation for He (P_He * V_ull = m_He * R_He * T_ull) yields:

m_He = P_He * V_ull / (R_He * T_ull)

dm_He/dt = (1 / R_He) * [ P_He * dV_ull/dt / T_ull
                         + V_ull * dP_He/dt / T_ull
                         - P_He * V_ull * dT_ull/dt / T_ull^2 ]

The terms dV_ull/dt, dP_He/dt, and dT_ull/dt can all be expressed analytically in terms of the current state and known quantities:

• dV_ull/dt = -dV_liq/dt = (m_dot_out - m_dot_in) / rho_liq (from liquid mass balance)
• dP_He/dt = -dP_N2/dt (since P_total is constant); dP_N2/dt is computed from the N2 ideal gas law with fixed m_N2_ull and known dV_ull/dt, dT_ull/dt
• dT_ull/dt is obtained from the ullage energy equation (which itself contains m_dot_He)

This creates a linear equation in m_dot_He that can be solved explicitly within each RHS evaluation. The key insight: dU_ull/dt = m_dot_He * h_in_He + Q_terms, and T_ull is a function of U_ull, so dT_ull/dt is linearly related to m_dot_He. Substituting into the dm_He/dt expression and solving for m_dot_He yields a closed-form formula with no finite differences, ensuring numerical stability with adaptive ODE solvers.

Implementation note: The analytical derivation will be implemented in tank_model.py as a dedicated method _solve_he_flow_rate() that returns m_dot_He as a function of the current state. This avoids the "algebraic loop" issue identified in review.

2.7 Heat Transfer Models

Liquid-to-ullage surface heat transfer:

Q_liq_to_ull = h_conv * A_surface * (T_liq - T_ull)

• A_surface: liquid surface area = cross-sectional area of cylinder = pi/4 * D^2
• h_conv: surface convective heat transfer coefficient [W/(m^2K)], configurable, default = 50 W/(m^2K)

Heat leak from environment (plugin interface): Total heat leak Q_leak is computed by the configured HeatLeakModel, then distributed to liquid and ullage zones by wetted area ratio:

Q_leak = heat_leak_model.compute(T_inner, T_env)
Q_leak_liq = Q_leak * A_wet / A_total
Q_leak_ull = Q_leak * A_dry / A_total

Where:

• A_wet = A_bottom + pi * D * liquid_level (bottom cap + wetted side wall)
• A_dry = A_top + pi * D * (H - liquid_level) (top cap + dry side wall)
• A_total = A_wet + A_dry
• T_inner passed to the model is a weighted average or conservative choice (e.g., T_liq for wet, T_ull for dry -- or simplified to just T_liq since most heat goes into the liquid)

Simplification: for the heat leak model input, use T_liq as T_inner since the liquid dominates thermal mass.

3. Heat Leak Interface

3.1 Base Class

class HeatLeakModel:
    def compute(self, T_inner: float, T_env: float) -> float:
        """Return total heat leak Q [W], positive = heat flows into tank."""
        raise NotImplementedError

3.2 MLI (Vacuum Multi-Layer Insulation)

class MLIHeatLeak(HeatLeakModel):
    def __init__(self, A_total, q_mli=1.0):
        """
        A_total: total tank surface area [m^2]
        q_mli: specific heat flux [W/m^2], default 1.0 W/m^2 (typical MLI performance)
        """

Computes: Q = A_total * q_mli

Note: MLI heat flux is largely independent of temperature difference in the typical operating range, so q_mli is a fixed parameter.

3.3 Foam/Wrap Insulation

class FoamHeatLeak(HeatLeakModel):
    def __init__(self, A_total, k_eff, delta):
        """
        A_total: total tank surface area [m^2]
        k_eff: effective thermal conductivity [W/(m*K)]
               Can be a float (constant) or a callable k_eff(T) -> float
               that returns conductivity as a function of temperature.
               At cryogenic temperatures, k varies significantly with T.
        delta: insulation thickness [m]
        """

Computes: Q = A_total * k(T_mean) * (T_env - T_inner) / delta

Where T_mean = (T_env + T_inner) / 2 when k_eff is a function, or simply uses the constant value when k_eff is a float.

4. Tank Geometry

Cylindrical tank:

• Total volume: V = 420.1 L = 0.4201 m^3
• Height: H = 0.5 m
• Cross-sectional area: A = V / H = 0.8402 m^2
• Diameter: D = sqrt(4*A/pi) = 1.034 m
• Side area: A_side = pi * D * H = 1.625 m^2
• Top area = Bottom area = A = 0.8402 m^2
• Total surface area: A_total = A_side + 2*A = 3.306 m^2

Liquid level at any time:

liquid_level = V_liq / A = (m_liq / rho_liq) / A
fill_fraction = liquid_level / H

5. Initial Conditions

Quantity	Value	Notes
T_liq(0)	78 K	Given
T_ull(0)	78 K	Given (same as liquid initially)
P_total	0.17 MPa	Constant throughout
Ullage fraction	30%	Gas pocket volume / total volume
V_liq(0)	0.2941 m^3	70% of 0.4201
V_ull(0)	0.1260 m^3	30% of 0.4201
m_liq(0)	rho_LN2(78K, 0.17MPa) * 0.2941	From CoolProp
U_liq(0)	m_liq(0) * u_LN2(78K, 0.17MPa)	Specific internal energy from CoolProp
m_N2_ull	rho_N2_vapor(78K, P_N2_sat) * V_ull(0)	N2 vapor at initial conditions, FIXED for all time
P_N2(0)	N2 saturation pressure at 78K	From CoolProp
P_He(0)	P_total - P_N2(0)	Helium makes up the pressure difference
m_He(0)	P_He(0) * V_ull(0) / (R_He * 78)	Ideal gas for He
U_ull(0)	m_N2_ull * u_N2_vapor(78K) + m_He(0) * u_He(78K)	Combined internal energy

6. Inlet/Outlet Conditions

Port	Flow rate	Temperature	Pressure	Fluid
LN2 inlet	1.144 kg/s	77 K	0.17 MPa	Liquid nitrogen
LN2 outlet	1.1895 kg/s	T_liq (computed)	0.17 MPa	Liquid nitrogen
He inlet	m_dot_He (computed)	100 K	0.17 MPa	Helium gas

Net liquid drain rate: 0.0455 kg/s.

7. Numerical Method

• Time integration: scipy.integrate.solve_ivp with RK45 (adaptive Runge-Kutta)
• Tolerances: rtol = 1e-8, atol = 1e-10
• Simulation duration: t_end = 3600 s (1 hour)
• Dense output: enabled for smooth interpolation of results
• State vector: y = [m_liq, U_liq, U_ull] (3 components)

RHS function evaluation at each call:

1. Unpack y -> (m_liq, U_liq, U_ull)
2. Compute T_liq from (m_liq, U_liq) via CoolProp
3. Compute V_liq, V_ull, liquid_level from geometry
4. Compute m_He from constant-pressure constraint
5. Compute T_ull from (U_ull, m_N2_ull, m_He) -- iterative or CoolProp
6. Compute all heat transfer terms (Q_liq_to_ull, Q_leak_liq, Q_leak_ull)
7. Compute enthalpy terms for inlets/outlet
8. Assemble and return dy/dt = [dm_liq/dt, dU_liq/dt, dU_ull/dt]

8. Code Structure

src/cryo_tank/
    __init__.py
    config.py          # Tank parameters, inlet/outlet conditions, simulation control
    properties.py      # CoolProp wrappers for N2 and He properties
    heat_leak.py       # HeatLeakModel base + MLIHeatLeak + FoamHeatLeak
    tank_model.py      # CryoTank class: geometry, state, rhs(), derived quantities
    solver.py          # run(tank, t_end) -> history dict
    output.py          # Plotting and reporting
    main.py            # Entry point

Module dependencies:

main.py -> config.py, tank_model.py, solver.py, output.py
tank_model.py -> properties.py, heat_leak.py
solver.py -> tank_model.py (calls tank.rhs)
output.py -> (only depends on history data dict)

8.1 CoolProp Performance Optimization (properties.py)

CoolProp Python calls can be slow when invoked thousands of times in an ODE RHS. The following optimizations are mandatory in properties.py:

1. Use CoolProp.AbstractState: Create persistent AbstractState objects for N2 and He at module load time. Reuse them across calls (avoid per-call overhead).

2. Lookup table with interpolation: For the dominant hot-path properties (LN2 density, enthalpy, internal energy at P = 0.17 MPa as a function of T), build a 1D interpolation table at startup covering the expected temperature range (e.g., 70-90 K for liquid, 70-300 K for gas). Use numpy.interp for fast evaluation. Fall back to CoolProp only when T is outside the table range.

3. He as ideal gas: Helium at 0.17 MPa and 78-300 K is well-described by the ideal gas law. Use analytical expressions (cp, cv, h, u) instead of CoolProp for He wherever possible to avoid unnecessary library calls.

9. Output Quantities

All recorded as time series:

Output	Symbol	Unit
Time	t	s
Liquid temperature	T_liq	K
Ullage temperature	T_ull	K
Liquid mass	m_liq	kg
Ullage N2 vapor mass	m_N2_ull	kg (constant)
Helium mass	m_He	kg
Helium flow rate	m_dot_He	kg/s
Liquid volume / level	V_liq, liquid_level	m^3, m
Fill fraction	fill_fraction	--
LN2 outlet temperature	T_out = T_liq	K
Heat leak total	Q_leak	W
Heat leak to liquid	Q_leak_liq	W
Heat leak to ullage	Q_leak_ull	W
Surface heat transfer	Q_liq_to_ull	W
Tank pressure	P_total	Pa (constant 0.17 MPa)

Output files:

• Time-series plots (PNG): T_liq & T_ull vs t, liquid level vs t, m_dot_He vs t, heat fluxes vs t
• History data (NPZ): all time series for post-processing
• Summary report (HTML): key parameters and embedded figures

10. Constraints, Limits, and Robustness

10.1 Normal Operating Limits

• Tank must not be overpressurized: P_total <= 0.8 MPa (max bearing pressure). Should raise warning/error if pressure constraint cannot be maintained.
• Liquid level must remain >= 0. When m_liq reaches 0, simulation should stop (use solve_ivp events mechanism to detect zero-crossing).

10.2 Near-Full Tank (fill_fraction > 95%)

When the ullage volume becomes very small, pressure sensitivity to mass/temperature changes grows dramatically. This can cause ODE solver instability.

Handling: when fill_fraction > 0.95, log a warning. The solver should still function because the analytical m_dot_He derivation avoids the finite-difference instability. If the solver fails to converge, reduce rtol/atol or switch to a stiffer solver (e.g., Radau).

10.3 Helium Backflow (m_dot_He < 0)

If the constant-pressure constraint computes m_dot_He < 0 (meaning pressure is too high and He should flow out), this indicates the operating regime has changed (e.g., heat leak is raising ullage temperature/pressure faster than liquid draining creates space). Two modes:

• Strict mode (default): Clamp m_dot_He = 0, allow pressure to drift above P_target. Log a warning with the overpressure magnitude. If P > 0.8 MPa (max bearing pressure), terminate simulation with an error.
• Vent mode (future): Add a pressure relief mechanism. Not implemented in v1.

10.4 Model Applicability (No Mass Transfer Assumption)

The current model assumes no evaporation/condensation between liquid and gas zones. This is valid when:

• Liquid temperature remains well below the saturation temperature at 0.17 MPa (~83.7 K)
• The net liquid drain is fast relative to temperature rise from heat leak

If T_liq approaches saturation temperature, the model will log a warning: "T_liq approaching saturation (83.7 K); evaporation effects may be significant."

Future extension: add Hertz-Knudsen evaporation model as an optional feature.

11. Testing Strategy

• Unit tests for properties.py: verify CoolProp wrappers return physically reasonable values
• Unit tests for heat_leak.py: verify MLI and Foam models compute correct Q for known inputs
• Unit tests for tank_model.py: verify geometry calculations, initial state, RHS evaluation
• Integration test: short run (10s), verify mass conservation (liquid mass change = net flow * dt)
• Steady-state test: with zero net flow and zero heat leak, verify state remains constant