402 lines
16 KiB
Markdown
402 lines
16 KiB
Markdown
# 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
|