# 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^2*K)], configurable, default = 50 W/(m^2*K)

**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

```python
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)

```python
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

```python
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