commit 0f57ba94f39c47f6d3bce37ba28e659fc5c297f8 Author: ljz <425868052@qq.com> Date: Wed Jun 3 15:41:04 2026 +0800 代码仓库移植 diff --git a/.codex b/.codex new file mode 100644 index 0000000..e69de29 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6fde6bb --- /dev/null +++ b/.gitignore @@ -0,0 +1,31 @@ +# OS +.DS_Store +Thumbs.db + +# Editors +.vscode/ +.idea/ +*.swp +*.swo + +# Python +__pycache__/ +*.pyc +.venv/ +venv/ + +# Build +build/ +dist/ +*.o +*.mod +*.out + +# Logs and temp +*.log +tmp/ +.cache/ + +# Results +results/* +!results/.gitkeep diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..901e90b --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,23 @@ +# Repository Guidelines + +## Project Structure & Module Organization +Core simulation code lives in `src/`. The base pipe blowdown model is split across `src/tank.py`, `src/pipe.py`, `src/riemann.py`, `src/solver.py`, and `src/output.py`, with `src/main.py` as the entry point. The cryogenic tank variant is isolated under `src/cryo_tank/`. Tests mirror that layout in `tests/` and `tests/cryo_tank/`. Use `cases/` for scenario inputs, `docs/` for design notes and generated reports, `scripts/` for utilities, and `results/` for run artifacts. + +## Build, Test, and Development Commands +There is no checked-in build system; run modules directly from the repository root. + +- `python3 src/main.py`: run the 0D-1D tank-pipe simulation and write plots/reports to `results/`. +- `python3 src/cryo_tank/main.py`: run the cryogenic LN2 tank simulation. +- `pytest -q`: run the full test suite. +- `pytest -q tests/test_integration.py`: run the main conservation tests only. +- `pytest -q tests/cryo_tank`: run the cryogenic tank test subset. +- `python3 generate_doc.py`: regenerate the Word technical document in `docs/` (requires `python-docx`). + +## Coding Style & Naming Conventions +Follow the existing Python style: 4-space indentation, `snake_case` for functions/modules, `PascalCase` for classes, and `UPPER_SNAKE_CASE` for configuration constants. Keep modules focused on one responsibility and prefer small helper functions like `_total_mass`. No formatter or linter config is checked in, so match current PEP 8-oriented style and keep imports simple and explicit. + +## Testing Guidelines +Tests use `pytest`, with `conftest.py` adding `src/` to `PYTHONPATH`. Name files `test_*.py` and keep related scenarios grouped by subsystem, for example `tests/test_pipe.py` or `tests/cryo_tank/test_integration.py`. Preserve the current emphasis on physical invariants: conservation, steady-state behavior, and plausible output trends. Add targeted regression tests whenever solver logic, property models, or boundary flux handling changes. + +## Commit & Pull Request Guidelines +Recent history follows Conventional Commit style: `feat: ...`, `feat(cryo_tank): ...`, `test(cryo_tank): ...`, `docs: ...`. Keep subjects imperative and concise. PRs should state the simulated scenario affected, summarize numerical or API changes, list test commands run, and attach updated plots or report outputs when behavior or post-processing changes. Avoid committing large generated files unless they are the point of the change. diff --git a/README.md b/README.md new file mode 100644 index 0000000..ec826c5 --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +# pipe-system-simulation-test + +Test repository for pipe system simulation work. + +## Structure + +- `docs/` design notes and reports +- `src/` source code +- `cases/` test cases and input data +- `scripts/` utility scripts +- `results/` generated outputs and post-processing summaries + +## Status + +Initial repository scaffold created. diff --git a/cases/README.md b/cases/README.md new file mode 100644 index 0000000..789d201 --- /dev/null +++ b/cases/README.md @@ -0,0 +1,3 @@ +# cases + +Simulation cases, configuration files, and input datasets. diff --git a/conftest.py b/conftest.py new file mode 100644 index 0000000..e7097f5 --- /dev/null +++ b/conftest.py @@ -0,0 +1,6 @@ +# conftest.py +import os +import sys + +_HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, os.path.join(_HERE, "src")) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..db65c75 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,3 @@ +# docs + +Project notes, design documents, and reports. diff --git a/docs/pipe_system_simulation_technical_doc.docx b/docs/pipe_system_simulation_technical_doc.docx new file mode 100644 index 0000000..14ad037 Binary files /dev/null and b/docs/pipe_system_simulation_technical_doc.docx differ diff --git a/docs/superpowers/plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md b/docs/superpowers/plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md new file mode 100644 index 0000000..b485e67 --- /dev/null +++ b/docs/superpowers/plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md @@ -0,0 +1,1426 @@ +# 0D–1D Tank-Pipe Blowdown MVP 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 simulation of a high-pressure tank blowing down through a 1 m pipe (D=5 mm) into a low-pressure tank, using a 0D lumped-parameter model for each tank coupled to a 1D finite-volume Euler solver on the pipe via ghost-cell + HLL Riemann fluxes. + +**Architecture:** 7 Python source modules (`config`, `riemann`, `tank`, `pipe`, `solver`, `output`, `main`) under `src/`, with `tests/` for pytest-based unit + integration tests. Test-driven workflow: each module is built by writing failing tests first, then the minimal implementation that passes them, then committing. The solver enforces mass/energy conservation via "flux doubling" — the same boundary flux vector updates both the pipe boundary cell and the adjacent tank. + +**Tech Stack:** Python 3, NumPy, Matplotlib (with `PillowWriter` for GIF animation, avoiding ffmpeg dependency), pytest. + +**Spec reference:** `docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md` + +--- + +## File Layout (Target State) + +``` +pipe-system-simulation-test/ +├── src/ +│ ├── config.py # constants (gas props, geometry, ICs, numerics) +│ ├── riemann.py # hll_flux(W_L, W_R, gamma) pure function +│ ├── tank.py # Tank class (0D, state = mass + total internal energy U) +│ ├── pipe.py # Pipe class (1D, conservative W field, step()) +│ ├── solver.py # run(tank1, tank2, pipe, t_end, cfl, ...) time loop +│ ├── output.py # save_history / plot_* / make_pipe_animation +│ └── main.py # entry: assemble → run → sanity check → save +├── tests/ +│ ├── test_riemann.py # 3 tests +│ ├── test_tank.py # 3 tests +│ ├── test_pipe.py # 2 tests +│ └── test_integration.py # 2 tests (total mass / energy conservation) +├── conftest.py # sys.path hack so tests/ can import from src/ +├── results/ # runtime outputs (.npz, .png, .gif) — created at runtime +└── docs/superpowers/ + ├── specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md + └── plans/2026-04-14-0d-1d-tank-pipe-blowdown-mvp.md (this file) +``` + +The `src/README.md`, `tests/`, `results/README.md`, `cases/README.md`, `scripts/README.md` existing scaffolding files are untouched; only `src/README.md` and friends are left as-is. + +--- + +## Pre-flight: Dependencies + +Verify Python 3 + required packages before starting any task. **Run this once at the beginning.** + +```bash +python3 -c "import numpy, matplotlib, pytest; print(numpy.__version__, matplotlib.__version__, pytest.__version__)" +``` + +Expected: three version strings printed with no import error. + +If `matplotlib` or `pytest` is missing: +```bash +pip install --user numpy matplotlib pytest +``` + +--- + +## Task 1: Project Skeleton + +Create the directory structure, `conftest.py` path shim, and `results/` directory. No Python code yet — just scaffolding so subsequent tasks can import cleanly. + +**Files:** +- Create: `conftest.py` +- Create: `tests/` (directory only) +- Create: `results/` (directory only; add a `.gitkeep` so git tracks the folder) + +- [ ] **Step 1: Create `conftest.py` at project root** + +Write this file exactly (path shim so `tests/test_*.py` can do `from riemann import ...`): + +```python +# conftest.py +import os +import sys + +_HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, os.path.join(_HERE, "src")) +``` + +- [ ] **Step 2: Create empty `tests/` directory** + +```bash +mkdir -p tests +``` + +- [ ] **Step 3: Create `results/` directory with `.gitkeep`** + +```bash +mkdir -p results +touch results/.gitkeep +``` + +- [ ] **Step 4: Verify pytest can discover the (empty) tests directory** + +```bash +pytest tests/ -v +``` + +Expected output: `no tests ran in X.XXs` (exit code 5 is fine — it means "no tests collected", which is what we want at this stage). + +- [ ] **Step 5: Commit** + +```bash +git add conftest.py tests results/.gitkeep +git commit -m "scaffold: add conftest, tests/, results/ for MVP" +``` + +--- + +## Task 2: `src/config.py` — Constants & Scenario + +A pure data module with all physical/numerical constants. Includes `assert` statements for parameter validation (per spec §6.1). No functions, no classes, no side effects. + +**Files:** +- Create: `src/config.py` + +- [ ] **Step 1: Write `src/config.py`** + +```python +# src/config.py +""" +Physical, geometric, and numerical constants for the 0D-1D tank-pipe +blowdown MVP. Pure data module — no functions, no side effects. +""" + +# ---------- Gas properties (ideal air-like) ---------- +GAMMA = 1.4 +R_GAS = 287.0 # J / (kg K) + +# ---------- High-pressure tank (upstream, Tank 1) ---------- +V1 = 5.0 # m^3 +P1_INIT = 10e6 # Pa (10 MPa) +T1_INIT = 300.0 # K + +# ---------- Low-pressure tank (downstream, Tank 2) ---------- +V2 = 10.0 # m^3 +P2_INIT = 101325.0 # Pa (1 atm) +T2_INIT = 300.0 # K + +# ---------- Pipe geometry ---------- +L = 1.0 # m +D = 5e-3 # m (5 mm) +N_CELLS = 20 # number of finite-volume cells + +# ---------- Simulation control ---------- +T_END = 0.1 # s +CFL = 0.5 + +# ---------- Output & animation ---------- +ANIMATION_STRIDE = 10 # keep every Nth frame in the GIF +OUTPUT_DIR = "results" + +# ---------- Parameter validation (per spec §6.1) ---------- +assert GAMMA > 1, "GAMMA must be > 1" +assert R_GAS > 0, "R_GAS must be > 0" +assert V1 > 0 and V2 > 0, "tank volumes must be > 0" +assert L > 0, "L must be > 0" +assert D > 0, "D must be > 0" +assert N_CELLS >= 2, "N_CELLS must be >= 2" +assert P1_INIT > 0 and P2_INIT > 0, "initial pressures must be > 0" +assert T1_INIT > 0 and T2_INIT > 0, "initial temperatures must be > 0" +assert 0 < CFL <= 1, "CFL must be in (0, 1]" +assert T_END > 0, "T_END must be > 0" +assert ANIMATION_STRIDE >= 1, "ANIMATION_STRIDE must be >= 1" +``` + +- [ ] **Step 2: Verify the module imports cleanly and assertions pass** + +```bash +python3 -c "import sys; sys.path.insert(0,'src'); import config; print(config.GAMMA, config.N_CELLS, config.T_END)" +``` + +Expected output: `1.4 20 0.1` + +- [ ] **Step 3: Commit** + +```bash +git add src/config.py +git commit -m "feat(config): add physical and numerical constants" +``` + +--- + +## Task 3: `src/riemann.py` — HLL Flux (TDD) + +HLL Riemann solver as a pure function. Written TDD: 3 tests first, then implementation. + +**Files:** +- Create: `tests/test_riemann.py` +- Create: `src/riemann.py` + +- [ ] **Step 1: Write the failing tests** + +Create `tests/test_riemann.py`: + +```python +# tests/test_riemann.py +import numpy as np +from riemann import hll_flux + + +GAMMA = 1.4 + + +def _to_conservative(rho, u, P, gamma): + """(rho, u, P) -> [rho, rho*u, rho*E] where E = e + u^2/2.""" + return np.array([ + rho, + rho * u, + P / (gamma - 1) + 0.5 * rho * u ** 2, + ]) + + +def _physical_flux(W, gamma): + """F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)].""" + rho = W[0] + u = W[1] / rho + E = W[2] + P = (gamma - 1) * (E - 0.5 * rho * u ** 2) + return np.array([rho * u, rho * u ** 2 + P, u * (E + P)]) + + +def test_hll_identical_states_returns_physical_flux(): + """ + When W_L == W_R, HLL must return the exact physical flux F(W) + with zero numerical dissipation (the (W_R - W_L) term vanishes). + """ + W = _to_conservative(rho=1.2, u=50.0, P=2.5e5, gamma=GAMMA) + F = hll_flux(W, W, GAMMA) + expected = _physical_flux(W, GAMMA) + assert np.allclose(F, expected, rtol=1e-12), f"F={F}, expected={expected}" + + +def test_hll_equal_pressure_equal_energy_gives_exact_pressure_flux(): + """ + Two stationary states (u=0) with same P but different rho: + - Both have the same energy density E = P/(gamma-1), so the HLL + (W_R - W_L)[2] term vanishes -> exact zero energy flux. + - F_L[1] = F_R[1] = P, and W_L[1] = W_R[1] = 0, so the momentum + flux is exactly P. + - The mass flux is NOT exactly zero for HLL (the density jump + triggers the (W_R - W_L)[0] dissipation term) — this is a known + HLL limitation for stationary contact discontinuities. We do + not assert on F[0] here. + """ + W_L = _to_conservative(rho=10.0, u=0.0, P=1e5, gamma=GAMMA) + W_R = _to_conservative(rho=1.0, u=0.0, P=1e5, gamma=GAMMA) + F = hll_flux(W_L, W_R, GAMMA) + assert abs(F[1] - 1e5) < 1e-6, f"momentum flux should equal P=1e5, got {F[1]}" + assert abs(F[2]) < 1e-8, f"energy flux should be exactly 0, got {F[2]}" + + +def test_hll_sod_shock_tube_directional_fluxes_all_positive(): + """ + Classical Sod initial values: + left: (rho, u, P) = (1.0, 0, 1.0) + right: (rho, u, P) = (0.125, 0, 0.1) + The pressure gradient drives flow from left to right, so the + HLL flux at the interface should have all three components + strictly positive: + F[0] > 0 : mass flux rightward + F[1] > 0 : momentum flux rightward + F[2] > 0 : energy flux rightward + """ + W_L = _to_conservative(rho=1.0, u=0.0, P=1.0, gamma=GAMMA) + W_R = _to_conservative(rho=0.125, u=0.0, P=0.1, gamma=GAMMA) + F = hll_flux(W_L, W_R, GAMMA) + assert F[0] > 0, f"expected positive mass flux, got {F[0]}" + assert F[1] > 0, f"expected positive momentum flux, got {F[1]}" + assert F[2] > 0, f"expected positive energy flux, got {F[2]}" +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +pytest tests/test_riemann.py -v +``` + +Expected: 3 failures with `ModuleNotFoundError: No module named 'riemann'`. + +- [ ] **Step 3: Write the minimal implementation** + +Create `src/riemann.py`: + +```python +# src/riemann.py +""" +HLL Riemann flux for the 1D compressible Euler equations. + +Conservative variable vector: W = [rho, rho*u, rho*E] + where E = e + u^2/2 is specific total energy, + e = P / (rho * (gamma - 1)) is specific internal energy. + +Physical flux: F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)] +""" + +import numpy as np + + +def hll_flux(W_L, W_R, gamma): + """ + Compute the HLL numerical flux at the interface between two states. + + Parameters + ---------- + W_L, W_R : array-like of shape (3,) + Left and right conservative state vectors. + gamma : float + Ratio of specific heats. + + Returns + ------- + np.ndarray of shape (3,) + HLL numerical flux vector. + + Raises + ------ + ValueError + If either state has non-positive density or pressure. + """ + # --- Recover primitives from left state --- + rho_L = W_L[0] + if rho_L <= 0: + raise ValueError(f"Non-positive density in W_L: rho={rho_L}, W_L={W_L}") + u_L = W_L[1] / rho_L + E_L = W_L[2] # total energy density (= rho * E) + p_L = (gamma - 1) * (E_L - 0.5 * rho_L * u_L ** 2) + if p_L <= 0: + raise ValueError(f"Non-positive pressure in W_L: p={p_L}, W_L={W_L}") + a_L = np.sqrt(gamma * p_L / rho_L) + + # --- Recover primitives from right state --- + rho_R = W_R[0] + if rho_R <= 0: + raise ValueError(f"Non-positive density in W_R: rho={rho_R}, W_R={W_R}") + u_R = W_R[1] / rho_R + E_R = W_R[2] + p_R = (gamma - 1) * (E_R - 0.5 * rho_R * u_R ** 2) + if p_R <= 0: + raise ValueError(f"Non-positive pressure in W_R: p={p_R}, W_R={W_R}") + a_R = np.sqrt(gamma * p_R / rho_R) + + # --- Physical fluxes at left/right states --- + F_L = np.array([ + rho_L * u_L, + rho_L * u_L ** 2 + p_L, + u_L * (E_L + p_L), + ]) + F_R = np.array([ + rho_R * u_R, + rho_R * u_R ** 2 + p_R, + u_R * (E_R + p_R), + ]) + + # --- HLL wave-speed estimates --- + S_L = min(u_L - a_L, u_R - a_R) + S_R = max(u_L + a_L, u_R + a_R) + + # --- HLL flux, piecewise on wave configuration --- + if S_L >= 0: + return F_L + if S_R <= 0: + return F_R + W_L_arr = np.asarray(W_L, dtype=float) + W_R_arr = np.asarray(W_R, dtype=float) + return (S_R * F_L - S_L * F_R + S_L * S_R * (W_R_arr - W_L_arr)) / (S_R - S_L) +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_riemann.py -v +``` + +Expected: `3 passed`. + +- [ ] **Step 5: Commit** + +```bash +git add src/riemann.py tests/test_riemann.py +git commit -m "feat(riemann): HLL flux with 3 unit tests" +``` + +--- + +## Task 4: `src/tank.py` — 0D Tank (TDD) + +Tank class with `(mass, U)` as primary state and `(rho, T, P)` as derived properties. Written TDD: 3 tests first, then implementation. + +**Files:** +- Create: `tests/test_tank.py` +- Create: `src/tank.py` + +- [ ] **Step 1: Write the failing tests** + +Create `tests/test_tank.py`: + +```python +# tests/test_tank.py +import numpy as np +import pytest +from tank import Tank + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def test_tank_initial_state_matches_ideal_gas(): + """ + Given P, T, V construct a Tank; mass and U should match ideal gas: + rho = P / (R * T) + mass = rho * V + U = P * V / (gamma - 1) (since u=0 inside tank) + T = (U / mass) * (gamma - 1) / R (round-trip) + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + rho_expected = 10e6 / (R_GAS * 300.0) + assert abs(tank.rho - rho_expected) < 1e-9 + assert abs(tank.mass - rho_expected * 5.0) < 1e-6 + assert abs(tank.U - 10e6 * 5.0 / (GAMMA - 1)) < 1e-3 + assert abs(tank.P - 10e6) < 1e-3 + assert abs(tank.T - 300.0) < 1e-9 + + +def test_tank_ghost_state_is_stagnation_conservative_vector(): + """ + ghost_state() should return [rho, 0, P/(gamma-1)] + (u_ghost = 0 per spec §2.3, so total energy density equals + internal energy density = P/(gamma-1)). + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + g = tank.ghost_state() + assert g.shape == (3,) + assert abs(g[0] - tank.rho) < 1e-12 + assert g[1] == 0.0 + assert abs(g[2] - 10e6 / (GAMMA - 1)) < 1e-3 + + +def test_tank_apply_flux_outflow_reduces_mass_energy_and_pressure(): + """ + apply_flux(mdot, edot, dt, sign=-1) should subtract mdot*dt from mass + and edot*dt from U. Derived P should decrease correspondingly. + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + mass_before = tank.mass + U_before = tank.U + P_before = tank.P + mdot = 1.0 # kg/s + edot = 5e5 # J/s (enthalpy rate) + dt = 1e-3 + tank.apply_flux(mdot=mdot, edot=edot, dt=dt, sign=-1) + assert abs(tank.mass - (mass_before - mdot * dt)) < 1e-12 + assert abs(tank.U - (U_before - edot * dt)) < 1e-9 + assert tank.P < P_before +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +pytest tests/test_tank.py -v +``` + +Expected: 3 failures with `ModuleNotFoundError: No module named 'tank'`. + +- [ ] **Step 3: Write the implementation** + +Create `src/tank.py`: + +```python +# src/tank.py +""" +0D lumped-parameter tank for ideal gas. The tank's *primary* state is +(mass, U) where U is total internal energy in joules. Pressure, temperature, +and density are derived properties computed from (mass, U) on demand, so +they are always consistent with the conservation-law updates. + +Conservation laws (u=0 inside tank): + dm/dt = mdot_in (mass) + dU/dt = Hdot_in = mdot_in * h_t,in (energy, open-system first law) + +where h_t is specific total enthalpy. When the tank couples to a 1D pipe +through the HLL boundary flux, flux[0]*A = mdot and flux[2]*A = Hdot +automatically — see solver.py. +""" + +import numpy as np + + +class Tank: + def __init__(self, V, P_init, T_init, gamma, R_gas): + self.V = V + self.gamma = gamma + self.R = R_gas + rho = P_init / (R_gas * T_init) + self.mass = rho * V + # For u=0, total internal energy equals rho*e*V = P*V / (gamma-1) + self.U = P_init * V / (gamma - 1) + + @property + def rho(self): + return self.mass / self.V + + @property + def T(self): + return (self.U / self.mass) * (self.gamma - 1) / self.R + + @property + def P(self): + return self.rho * self.R * self.T + + def ghost_state(self): + """ + Return the conservative variable vector [rho, rho*u, rho*E] that + represents this tank as a ghost cell for the 1D pipe solver. + Since u_ghost = 0, rho*u = 0 and rho*E = P/(gamma-1). + """ + return np.array([self.rho, 0.0, self.P / (self.gamma - 1)]) + + def apply_flux(self, mdot, edot, dt, sign): + """ + Update (mass, U) from one time step of boundary flux. + + Parameters + ---------- + mdot : float + Mass flux across the interface in kg/s (already multiplied + by pipe cross-sectional area). Sign is the "outward normal" + convention of the pipe: positive = pipe-rightward. + edot : float + Total enthalpy rate in W (= flux[2] * A), same convention. + dt : float + Time-step size in seconds. + sign : int (+1 or -1) + Orientation for this tank. For an upstream tank whose gas + flows "out to the right" into the pipe, the HLL left-boundary + flux has mdot > 0, so sign = -1 (tank loses mass). + For a downstream tank receiving gas from the right boundary + with mdot > 0 entering, sign = +1. + + Raises + ------ + RuntimeError + If the tank's mass becomes non-positive after the update. + """ + self.mass += sign * mdot * dt + self.U += sign * edot * dt + if self.mass <= 0: + raise RuntimeError( + f"Tank mass non-positive after apply_flux: mass={self.mass}, " + f"mdot={mdot}, edot={edot}, dt={dt}, sign={sign}" + ) +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_tank.py -v +``` + +Expected: `3 passed`. + +- [ ] **Step 5: Commit** + +```bash +git add src/tank.py tests/test_tank.py +git commit -m "feat(tank): 0D tank with (mass, U) primary state and 3 tests" +``` + +--- + +## Task 5: `src/pipe.py` — 1D Pipe (TDD) + +Pipe class storing the conservative-variable field W of shape (3, N). `step()` takes left/right boundary fluxes and advances the field by one explicit Euler step. + +**Files:** +- Create: `tests/test_pipe.py` +- Create: `src/pipe.py` + +- [ ] **Step 1: Write the failing tests** + +Create `tests/test_pipe.py`: + +```python +# tests/test_pipe.py +import numpy as np +import pytest +from pipe import Pipe + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def test_pipe_uniform_initialization_all_cells_identical(): + """ + Initial state: uniform P, T, u=0 across all N cells. + All cells should have identical W. primitives() should round-trip + back to u=0 and P=P_init. + """ + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=1e5, T_init=300.0, + gamma=GAMMA, R_gas=R_GAS) + + # All cells identical + for i in range(1, pipe.N): + assert np.allclose(pipe.W[:, i], pipe.W[:, 0], rtol=1e-14) + + rho, u, P, a = pipe.primitives() + rho_expected = 1e5 / (R_GAS * 300.0) + assert np.allclose(u, 0.0) + assert np.allclose(P, 1e5, rtol=1e-10) + assert np.allclose(rho, rho_expected, rtol=1e-10) + + # Sound speed a = sqrt(gamma * P / rho) + a_expected = np.sqrt(GAMMA * 1e5 / rho_expected) + assert np.allclose(a, a_expected, rtol=1e-10) + + # Geometric sanity + assert pipe.dx == 1.0 / 20 + assert pipe.x_centers.shape == (20,) + assert abs(pipe.x_centers[0] - 0.025) < 1e-14 + assert abs(pipe.x_centers[-1] - 0.975) < 1e-14 + + +def test_pipe_uniform_state_plus_matching_boundary_flux_is_static(): + """ + If all cells have identical W (so all internal HLL fluxes are + identical to the physical flux F(W) = [0, P, 0] for u=0 state), + and we pass in boundary fluxes equal to [0, P, 0] as well, then + every difference (flux[i+1] - flux[i]) is zero, so W must not + change after one step. Verify to machine precision. + """ + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=1e5, T_init=300.0, + gamma=GAMMA, R_gas=R_GAS) + + W_before = pipe.W.copy() + # Boundary flux matching the static interior: [rho*u, rho*u^2+P, u*(E+P)] + # with u=0 -> [0, P, 0] + flux_boundary = np.array([0.0, 1e5, 0.0]) + pipe.step(flux_boundary, flux_boundary, dt=1e-5) + assert np.allclose(pipe.W, W_before, atol=1e-6, rtol=1e-12) +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +pytest tests/test_pipe.py -v +``` + +Expected: 2 failures with `ModuleNotFoundError: No module named 'pipe'`. + +- [ ] **Step 3: Write the implementation** + +Create `src/pipe.py`: + +```python +# src/pipe.py +""" +1D finite-volume pipe for compressible Euler equations: + dW/dt + dF(W)/dx = 0 + W = [rho, rho*u, rho*E], F = [rho*u, rho*u^2 + P, u*(rho*E + P)] + +Discretization: + - N uniform cells, cell-averaged piecewise-constant reconstruction + - HLL numerical flux at all interior interfaces + - Boundary (tank-side) interface fluxes are provided by the caller + via step(flux_L, flux_R, dt) +""" + +import numpy as np +from riemann import hll_flux + + +class Pipe: + def __init__(self, L, D, N, P_init, T_init, gamma, R_gas): + self.L = L + self.D = D + self.N = N + self.dx = L / N + self.area = np.pi * (D / 2) ** 2 + self.gamma = gamma + self.R = R_gas + self.x_centers = np.linspace(self.dx / 2, L - self.dx / 2, N) + + # Uniform initial state, u = 0 + rho = P_init / (R_gas * T_init) + E_density = P_init / (gamma - 1) # since u=0, total energy density = internal + + self.W = np.zeros((3, N)) + self.W[0, :] = rho + self.W[1, :] = 0.0 + self.W[2, :] = E_density + + def primitives(self): + """ + Return (rho, u, P, a) each of shape (N,), computed from W. + """ + rho = self.W[0, :] + u = self.W[1, :] / rho + P = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho * u ** 2) + a = np.sqrt(self.gamma * P / rho) + return rho, u, P, a + + def max_wave_speed(self): + """ + Return max over cells of |u| + a, used for CFL dt calculation. + """ + _, u, _, a = self.primitives() + return float(np.max(np.abs(u) + a)) + + def step(self, flux_L, flux_R, dt): + """ + Advance W by one explicit Euler step. The caller provides the + two boundary interface fluxes (with ghost states already folded + in); internal interface fluxes are computed here with HLL. + + Parameters + ---------- + flux_L, flux_R : np.ndarray of shape (3,) + Numerical fluxes at the leftmost and rightmost interfaces + (cell -1/2 and cell N-1/2, i.e. the tank-facing boundaries). + dt : float + Time-step size. + + Raises + ------ + RuntimeError + If the updated state has any non-positive density or pressure. + """ + N = self.N + W_snap = self.W.copy() + + # Internal fluxes: flux_int[:, k] is the flux at the interface + # between cell k and cell k+1, for k = 0 .. N-2 (total N-1 of them) + flux_int = np.zeros((3, N - 1)) + for k in range(N - 1): + flux_int[:, k] = hll_flux(W_snap[:, k], W_snap[:, k + 1], self.gamma) + + # First cell: left face = flux_L, right face = flux_int[:, 0] + self.W[:, 0] = W_snap[:, 0] - (dt / self.dx) * (flux_int[:, 0] - flux_L) + + # Interior cells: left face = flux_int[:, i-1], right face = flux_int[:, i] + for i in range(1, N - 1): + self.W[:, i] = W_snap[:, i] - (dt / self.dx) * (flux_int[:, i] - flux_int[:, i - 1]) + + # Last cell: left face = flux_int[:, N-2], right face = flux_R + self.W[:, N - 1] = W_snap[:, N - 1] - (dt / self.dx) * (flux_R - flux_int[:, N - 2]) + + # Physical-state sanity check + rho_new = self.W[0, :] + if np.any(rho_new <= 0): + bad = np.where(rho_new <= 0)[0] + raise RuntimeError( + f"Non-positive density after pipe step at cells {bad.tolist()}: " + f"rho={rho_new[bad].tolist()}" + ) + u_new = self.W[1, :] / rho_new + P_new = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho_new * u_new ** 2) + if np.any(P_new <= 0): + bad = np.where(P_new <= 0)[0] + raise RuntimeError( + f"Non-positive pressure after pipe step at cells {bad.tolist()}: " + f"P={P_new[bad].tolist()}" + ) +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_pipe.py -v +``` + +Expected: `2 passed`. + +- [ ] **Step 5: Commit** + +```bash +git add src/pipe.py tests/test_pipe.py +git commit -m "feat(pipe): 1D FV pipe with HLL internal fluxes and 2 tests" +``` + +--- + +## Task 6: `src/solver.py` + Integration Tests (TDD) + +The time-loop driver that orchestrates tank ↔ pipe coupling. The integration tests verify total mass and total energy are conserved to 1e-10 over a short run — the end-to-end correctness check for the "flux doubling" mechanism. + +**Files:** +- Create: `tests/test_integration.py` +- Create: `src/solver.py` + +- [ ] **Step 1: Write the failing integration tests** + +Create `tests/test_integration.py`: + +```python +# tests/test_integration.py +""" +End-to-end conservation tests: after a short run, total system mass +and total system energy must be conserved to machine precision. These +are the load-bearing tests for the "flux doubling" coupling mechanism. +""" +import numpy as np +import pytest +from tank import Tank +from pipe import Pipe +from solver import run + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def _build_scenario(): + """Default blowdown scenario at smaller scale — same physics, short run.""" + tank1 = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + tank2 = Tank(V=10.0, P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + return tank1, tank2, pipe + + +def _total_mass(tank1, tank2, pipe): + return tank1.mass + tank2.mass + float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx)) + + +def _total_energy(tank1, tank2, pipe): + return tank1.U + tank2.U + float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx)) + + +def test_total_mass_conserved_over_short_run(): + tank1, tank2, pipe = _build_scenario() + m_init = _total_mass(tank1, tank2, pipe) + run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5) + m_final = _total_mass(tank1, tank2, pipe) + rel_err = abs(m_final - m_init) / m_init + assert rel_err < 1e-10, ( + f"Total mass not conserved: m_init={m_init:.6e}, " + f"m_final={m_final:.6e}, rel_err={rel_err:.2e}" + ) + + +def test_total_energy_conserved_over_short_run(): + tank1, tank2, pipe = _build_scenario() + U_init = _total_energy(tank1, tank2, pipe) + run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5) + U_final = _total_energy(tank1, tank2, pipe) + rel_err = abs(U_final - U_init) / U_init + assert rel_err < 1e-10, ( + f"Total energy not conserved: U_init={U_init:.6e}, " + f"U_final={U_final:.6e}, rel_err={rel_err:.2e}" + ) +``` + +- [ ] **Step 2: Run the tests to verify they fail** + +```bash +pytest tests/test_integration.py -v +``` + +Expected: 2 failures with `ModuleNotFoundError: No module named 'solver'`. + +- [ ] **Step 3: Write the solver implementation** + +Create `src/solver.py`: + +```python +# src/solver.py +""" +Time-loop driver for the 0D-1D coupled tank-pipe simulation. + +Per time step (per spec §3.1): + 1. Compute CFL-limited dt from pipe's max wave speed + 2. Freeze ghost states from current tank states + 3. Compute two boundary HLL fluxes (left and right) + 4. Advance pipe by one step using those two fluxes (pipe.step handles + the internal fluxes itself) + 5. Advance both tanks using the SAME two boundary fluxes * area + -> this "flux doubling" is the mechanism that makes system mass + and energy strictly conserved to machine precision + 6. Advance time + 7. Append snapshot to history +""" +import numpy as np +from riemann import hll_flux + + +def run(tank1, tank2, pipe, t_end, cfl, verbose=False, log_every=100): + """ + Run the coupled tank-pipe simulation from t=0 to t=t_end. + + Parameters + ---------- + tank1, tank2 : Tank + Upstream and downstream tanks. tank1 connects to pipe.W[:, 0], + tank2 connects to pipe.W[:, -1]. + pipe : Pipe + 1D pipe instance with initial state already set. + t_end : float + End time in seconds. + cfl : float + CFL number in (0, 1]. + verbose : bool, default False + If True, print step-progress info every `log_every` steps. + log_every : int, default 100 + Logging interval when verbose=True. + + Returns + ------- + dict + History with keys 't', 'P1', 'T1', 'P2', 'T2' (all 1D arrays + of shape (n_steps,)), and 'W_hist' of shape (n_steps, 3, N). + """ + history = { + 't': [], + 'P1': [], 'T1': [], + 'P2': [], 'T2': [], + 'W_hist': [], + } + + t = 0.0 + step = 0 + + while t < t_end: + # --- Phase 1: CFL time step --- + a_max = pipe.max_wave_speed() + dt = cfl * pipe.dx / a_max + dt = min(dt, t_end - t) + if dt < 1e-12: + raise RuntimeError( + f"dt degenerate at step {step}: dt={dt:.3e}, a_max={a_max:.3e}" + ) + + # --- Phase 2: freeze tank ghost states (snapshot for this step) --- + W_ghost_L = tank1.ghost_state() + W_ghost_R = tank2.ghost_state() + + # --- Phase 3: two boundary fluxes (solver-level) --- + flux_L = hll_flux(W_ghost_L, pipe.W[:, 0], pipe.gamma) + flux_R = hll_flux(pipe.W[:, -1], W_ghost_R, pipe.gamma) + + # --- Phase 4: advance pipe (internal fluxes handled inside) --- + pipe.step(flux_L, flux_R, dt) + + # --- Phase 5: advance tanks with the SAME boundary fluxes * area --- + fL_A = flux_L * pipe.area + fR_A = flux_R * pipe.area + # Left boundary flux is "rightward positive"; tank1 loses that mass + tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=dt, sign=-1) + # Right boundary flux is "rightward positive"; tank2 gains that mass + tank2.apply_flux(mdot=fR_A[0], edot=fR_A[2], dt=dt, sign=+1) + + # --- Phase 6: advance time --- + t += dt + step += 1 + + # --- Phase 7: record history --- + history['t'].append(t) + history['P1'].append(tank1.P) + history['T1'].append(tank1.T) + history['P2'].append(tank2.P) + history['T2'].append(tank2.T) + history['W_hist'].append(pipe.W.copy()) + + if verbose and step % log_every == 0: + _, u, _, _ = pipe.primitives() + print( + f"step={step:6d} t={t:.5f} dt={dt:.2e} " + f"P1={tank1.P/1e6:7.4f}MPa P2={tank2.P/1e6:7.4f}MPa " + f"max|u|={float(np.max(np.abs(u))):7.1f}m/s" + ) + + if step == 0: + raise RuntimeError("solver.run() exited without taking any step") + + # Convert lists to arrays for downstream consumers + history['t'] = np.asarray(history['t']) + history['P1'] = np.asarray(history['P1']) + history['T1'] = np.asarray(history['T1']) + history['P2'] = np.asarray(history['P2']) + history['T2'] = np.asarray(history['T2']) + history['W_hist'] = np.stack(history['W_hist']) # shape (n_steps, 3, N) + + return history +``` + +- [ ] **Step 4: Run integration tests to verify they pass** + +```bash +pytest tests/test_integration.py -v +``` + +Expected: `2 passed`. + +- [ ] **Step 5: Run the full test suite to verify nothing regressed** + +```bash +pytest tests/ -v +``` + +Expected: `10 passed` total (3 riemann + 3 tank + 2 pipe + 2 integration). + +- [ ] **Step 6: Commit** + +```bash +git add src/solver.py tests/test_integration.py +git commit -m "feat(solver): coupled time loop with mass/energy conservation tests" +``` + +--- + +## Task 7: `src/output.py` — Persistence, Plotting, Animation + +Utility module for `.npz` persistence, static time-series plots, and the pipe evolution GIF. No unit tests per spec §7.3 (visual outputs verified by eye). + +**Files:** +- Create: `src/output.py` + +- [ ] **Step 1: Write `src/output.py`** + +```python +# src/output.py +""" +Output helpers: persistence (.npz), static plots (.png), animation (.gif). +Uses matplotlib's Agg backend so it works in headless environments. +The PillowWriter is used for GIF output to avoid an ffmpeg dependency. +""" +import os + +import numpy as np +import matplotlib +matplotlib.use("Agg") # headless-safe; must be set before pyplot import +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation, PillowWriter + + +def save_history(history, pipe, path, gamma, R_gas): + """ + Persist the full simulation history + pipe geometry to a .npz file. + + Loadable later with: + d = np.load("results/history.npz") + rho = d['W_hist'][:, 0, :] + u = d['W_hist'][:, 1, :] / rho + P = (d['gamma'] - 1) * (d['W_hist'][:, 2, :] - 0.5 * rho * u**2) + """ + dirname = os.path.dirname(path) + if dirname: + os.makedirs(dirname, exist_ok=True) + np.savez_compressed( + path, + t=history['t'], + P1=history['P1'], T1=history['T1'], + P2=history['P2'], T2=history['T2'], + W_hist=history['W_hist'], + x=pipe.x_centers, + dx=pipe.dx, + area=pipe.area, + gamma=gamma, + R_gas=R_gas, + ) + + +def plot_tank_pressure(history, path): + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(history['t'], history['P1'] / 1e6, label="Tank 1 (high pressure)") + ax.plot(history['t'], history['P2'] / 1e6, label="Tank 2 (low pressure)") + ax.set_xlabel("Time [s]") + ax.set_ylabel("Pressure [MPa]") + ax.set_title("Tank pressures vs time") + ax.grid(True) + ax.legend() + fig.tight_layout() + fig.savefig(path, dpi=120) + plt.close(fig) + + +def plot_tank_temperature(history, path): + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(history['t'], history['T1'], label="Tank 1 (high pressure)") + ax.plot(history['t'], history['T2'], label="Tank 2 (low pressure)") + 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 make_pipe_animation(history, pipe, path, gamma, R_gas, stride=10): + """ + Render a GIF of the pipe's P(x), u(x), T(x) evolution over time. + Uses PillowWriter so no ffmpeg is needed. + """ + t = history['t'] + W_hist = history['W_hist'] # (n_steps, 3, N) + x = pipe.x_centers + n_steps = W_hist.shape[0] + + # Frame indices: every `stride`th snapshot, plus the final one + frames = list(range(0, n_steps, stride)) + if frames[-1] != n_steps - 1: + frames.append(n_steps - 1) + + # Precompute primitives for all frames in one vectorized pass + rho = W_hist[:, 0, :] + u = W_hist[:, 1, :] / rho + P = (gamma - 1) * (W_hist[:, 2, :] - 0.5 * rho * u ** 2) + T = P / (rho * R_gas) + + fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) + + # Pressure subplot + line_P, = axes[0].plot(x, P[0] / 1e6) + axes[0].set_ylabel("P [MPa]") + axes[0].set_ylim(P.min() / 1e6 * 0.95, P.max() / 1e6 * 1.05) + axes[0].grid(True) + + # Velocity subplot + line_u, = axes[1].plot(x, u[0]) + axes[1].set_ylabel("u [m/s]") + u_min, u_max = float(u.min()), float(u.max()) + pad = max(1.0, 0.05 * (u_max - u_min) if u_max > u_min else 1.0) + axes[1].set_ylim(u_min - pad, u_max + pad) + axes[1].grid(True) + + # Temperature subplot + line_T, = axes[2].plot(x, T[0]) + axes[2].set_ylabel("T [K]") + axes[2].set_xlabel("x [m]") + axes[2].set_ylim(T.min() * 0.95, T.max() * 1.05) + axes[2].grid(True) + + title = fig.suptitle("") + + def update(frame_idx): + line_P.set_ydata(P[frame_idx] / 1e6) + line_u.set_ydata(u[frame_idx]) + line_T.set_ydata(T[frame_idx]) + title.set_text( + f"t = {t[frame_idx]:.5f} s (step {frame_idx + 1}/{n_steps})" + ) + return line_P, line_u, line_T, title + + anim = FuncAnimation(fig, update, frames=frames, interval=50, blit=False) + + dirname = os.path.dirname(path) + if dirname: + os.makedirs(dirname, exist_ok=True) + anim.save(path, writer=PillowWriter(fps=20)) + plt.close(fig) +``` + +- [ ] **Step 2: Smoke-test that the module imports without errors** + +```bash +python3 -c "import sys; sys.path.insert(0,'src'); import output; print('output OK:', dir(output))" +``` + +Expected: `output OK: [...] 'make_pipe_animation', 'plot_tank_pressure', 'plot_tank_temperature', 'save_history', ...` + +- [ ] **Step 3: Commit** + +```bash +git add src/output.py +git commit -m "feat(output): npz persistence, timeseries plots, pipe GIF animation" +``` + +--- + +## Task 8: `src/main.py` — Entry Point + Sanity Check + +Wires everything together: build objects from `config`, call `solver.run`, assert conservation, write all outputs. + +**Files:** +- Create: `src/main.py` + +- [ ] **Step 1: Write `src/main.py`** + +```python +# src/main.py +""" +Entry point: assemble tanks + pipe from config constants, run the solver, +verify total mass/energy conservation, persist history, and generate +plots + animation. + +Run from project root: + python3 src/main.py +""" +import os +import sys + +# Ensure imports work when running from project root +_HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, _HERE) + +import numpy as np + +from config import ( + GAMMA, R_GAS, + V1, P1_INIT, T1_INIT, + V2, P2_INIT, T2_INIT, + L, D, N_CELLS, + T_END, CFL, + ANIMATION_STRIDE, OUTPUT_DIR, +) +from tank import Tank +from pipe import Pipe +from solver import run +from output import ( + save_history, + plot_tank_pressure, + plot_tank_temperature, + make_pipe_animation, +) + + +def _total_mass(tank1, tank2, pipe): + pipe_mass = float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx)) + return tank1.mass + tank2.mass + pipe_mass + + +def _total_energy(tank1, tank2, pipe): + pipe_energy = float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx)) + return tank1.U + tank2.U + pipe_energy + + +def main(): + os.makedirs(OUTPUT_DIR, exist_ok=True) + + # --- Assemble --- + tank1 = Tank(V=V1, P_init=P1_INIT, T_init=T1_INIT, gamma=GAMMA, R_gas=R_GAS) + tank2 = Tank(V=V2, P_init=P2_INIT, T_init=T2_INIT, gamma=GAMMA, R_gas=R_GAS) + pipe = Pipe(L=L, D=D, N=N_CELLS, P_init=P2_INIT, T_init=T2_INIT, + gamma=GAMMA, R_gas=R_GAS) + + m_init = _total_mass(tank1, tank2, pipe) + U_init = _total_energy(tank1, tank2, pipe) + print(f"Initial total mass: {m_init:.6e} kg") + print(f"Initial total energy: {U_init:.6e} J") + print(f"Initial P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa") + print(f"Pipe: L={L} m, D={D*1e3:.1f} mm, N={N_CELLS} cells, dx={pipe.dx*1e3:.1f} mm") + print(f"Running to t_end={T_END} s with CFL={CFL}...") + print() + + # --- Run --- + history = run(tank1, tank2, pipe, + t_end=T_END, cfl=CFL, + verbose=True, log_every=200) + + n_steps = len(history['t']) + print() + print(f"Simulation complete: {n_steps} steps") + + # --- Conservation sanity check (per spec §6.1, §8) --- + m_final = _total_mass(tank1, tank2, pipe) + U_final = _total_energy(tank1, tank2, pipe) + rel_err_m = abs(m_final - m_init) / m_init + rel_err_U = abs(U_final - U_init) / U_init + print(f"Final total mass: {m_final:.6e} kg (rel err = {rel_err_m:.2e})") + print(f"Final total energy: {U_final:.6e} J (rel err = {rel_err_U:.2e})") + print(f"Final P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa") + assert rel_err_m < 1e-10, f"Total mass not conserved: rel_err={rel_err_m:.2e}" + assert rel_err_U < 1e-10, f"Total energy not conserved: rel_err={rel_err_U:.2e}" + + # --- Persist + visualize --- + save_history(history, pipe, + os.path.join(OUTPUT_DIR, "history.npz"), + GAMMA, R_GAS) + plot_tank_pressure(history, + os.path.join(OUTPUT_DIR, "tank_pressure.png")) + plot_tank_temperature(history, + os.path.join(OUTPUT_DIR, "tank_temperature.png")) + make_pipe_animation(history, pipe, + os.path.join(OUTPUT_DIR, "pipe_animation.gif"), + GAMMA, R_GAS, stride=ANIMATION_STRIDE) + + print(f"Outputs written to {OUTPUT_DIR}/") + print(f" - history.npz") + print(f" - tank_pressure.png") + print(f" - tank_temperature.png") + print(f" - pipe_animation.gif") + + +if __name__ == "__main__": + main() +``` + +- [ ] **Step 2: Verify the module imports cleanly (no run)** + +```bash +python3 -c "import sys; sys.path.insert(0,'src'); import main; print('main module OK')" +``` + +Expected: `main module OK` (importing must not trigger `main()` because of `if __name__ == '__main__'`). + +- [ ] **Step 3: Commit** + +```bash +git add src/main.py +git commit -m "feat(main): entry point with conservation sanity check" +``` + +--- + +## Task 9: End-to-End Verification Run + +Run the full program and verify all acceptance criteria from spec §8. + +**Files:** +- No new files; verifies existing ones. + +- [ ] **Step 1: Run the full simulation** + +```bash +python3 src/main.py +``` + +Expected behavior: +- Prints "Initial total mass", "Initial total energy", initial pressures +- Prints verbose step logs every 200 steps (roughly 10–15 log lines total) +- Prints "Simulation complete: NNNN steps" where NNNN ≈ 2000–3000 +- Prints "Final total mass" and "rel err" < 1e-10 +- Prints "Final total energy" and "rel err" < 1e-10 +- Prints "Outputs written to results/" +- Exit code 0 +- Total wall time < 30 seconds + +- [ ] **Step 2: Verify all 4 output files exist and are non-empty** + +```bash +ls -la results/history.npz results/tank_pressure.png results/tank_temperature.png results/pipe_animation.gif +``` + +Expected: all four files present, each with non-zero size (`history.npz` ≈ 1-3 MB, PNGs tens of KB, GIF 5-15 MB). + +- [ ] **Step 3: Verify the .npz file is loadable and has correct shapes** + +```bash +python3 -c " +import numpy as np +d = np.load('results/history.npz') +print('keys:', sorted(d.files)) +print('t shape:', d['t'].shape) +print('W_hist shape:', d['W_hist'].shape) +print('P1 range:', float(d['P1'].min())/1e6, '->', float(d['P1'].max())/1e6, 'MPa') +print('P2 range:', float(d['P2'].min())/1e6, '->', float(d['P2'].max())/1e6, 'MPa') +# Physical monotonicity +assert d['P1'][0] > d['P1'][-1], 'P1 should decrease' +assert d['P2'][0] < d['P2'][-1], 'P2 should increase' +print('monotonicity OK') +" +``` + +Expected: +- `keys: ['P1', 'P2', 'R_gas', 'T1', 'T2', 'W_hist', 'area', 'dx', 'gamma', 't', 'x']` +- `t shape: (NNNN,)` where NNNN ≈ 2000–3000 +- `W_hist shape: (NNNN, 3, 20)` +- `P1 range`: starts near 10, decreases toward some intermediate value +- `P2 range`: starts near 0.1, increases toward some intermediate value +- `monotonicity OK` + +- [ ] **Step 4: Run the full pytest suite one more time** + +```bash +pytest tests/ -v +``` + +Expected: `10 passed`. + +- [ ] **Step 5: Commit a `.gitignore` entry for the bulky result files** + +Check whether `.gitignore` already excludes `results/`: + +```bash +cat .gitignore +``` + +If `results/` is not already ignored, update `.gitignore` to include it (but keep `.gitkeep`): + +Add these two lines to `.gitignore`: + +``` +results/* +!results/.gitkeep +``` + +Then: + +```bash +git add .gitignore +git commit -m "chore: ignore runtime outputs under results/" +``` + +If `results/` was already handled, skip this step with `git status` confirming no .gitignore changes are staged. + +- [ ] **Step 6: Final status check** + +```bash +git status +git log --oneline -12 +``` + +Expected: +- `git status` shows a clean working tree +- `git log` shows the sequence of feature commits: scaffold → config → riemann → tank → pipe → solver → output → main → (optionally) .gitignore + +--- + +## Acceptance Criteria Recap (from spec §8) + +After Task 9 is complete, verify each item is satisfied: + +1. [x] `python3 src/main.py` runs to completion without exceptions → Task 9 Step 1 +2. [x] `pytest tests/ -v` → 10 passed → Task 9 Step 4 +3. [x] main.py sanity check: total mass & total energy rel err < 1e-10 → Task 9 Step 1 output +4. [x] `P1(t)` monotone decreasing, `P2(t)` monotone increasing → Task 9 Step 3 +5. [x] `pipe_animation.gif` is a valid GIF openable by standard viewers → Task 9 Step 2 + manual inspection +6. [x] `history.npz` loadable via `np.load` with all fields accessible → Task 9 Step 3 +7. [x] Wall time < 30 seconds → Task 9 Step 1 diff --git a/docs/superpowers/plans/2026-04-16-cryo-tank-module.md b/docs/superpowers/plans/2026-04-16-cryo-tank-module.md new file mode 100644 index 0000000..1c548d7 --- /dev/null +++ b/docs/superpowers/plans/2026-04-16-cryo-tank-module.md @@ -0,0 +1,1432 @@ +# 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" +``` diff --git a/docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md b/docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md new file mode 100644 index 0000000..8f502bb --- /dev/null +++ b/docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md @@ -0,0 +1,486 @@ +# 0D–1D 气瓶–管路耦合瞬态放气仿真 MVP + +**日期**:2026-04-14 +**状态**:草案,待评审 +**作者**:He Yunqin(通过 Claude Code / superpowers:brainstorming) + +--- + +## 1. 目标与范围 + +### 1.1 目标 + +实现一个可运行的瞬态仿真程序,模拟以下场景: + +- **高压气瓶** V₁ = 5 m³,P₁ = 10 MPa,T₁ = 300 K +- **低压气瓶** V₂ = 10 m³,P₂ = 0.101325 MPa(1 atm),T₂ = 300 K +- **连接管路** L = 1 m,D = 5 mm,划分为 N = 20 个有限体积单元 +- **工质**:理想气体,γ = 1.4,R = 287 J/(kg·K) +- **初始条件**:管路与低压瓶同压同温(P₂, T₂),高压瓶独立 +- **仿真时长**:t_end = 0.1 s(捕获开启瞬间的激波与膨胀波过程) + +程序应输出: +1. 双瓶 P(t)、T(t) 时间曲线图 +2. 管内 P(x)、u(x)、T(x) 演化动画(GIF) +3. 完整时间序列数据文件(`.npz`),支持离线加载任意时刻任意空间点的状态 + +### 1.2 MVP 范围边界 + +**本 MVP 做**: + +- 0D 集中参数气瓶 + 1D 可压缩欧拉方程管路 +- 虚网格(Ghost Cell)+ HLL Riemann 求解器耦合 +- 一阶空间重构 + 一阶显式欧拉时间推进 +- CFL 动态步长控制 +- 基本 pytest 单元测试(共 10 条) + +**本 MVP 明确不做**: + +- 摩擦源项(Darcy–Weisbach 或 Fanning) +- 壁面传热(绝热假设) +- 二阶空间格式(MUSCL / limiter) +- 二阶时间格式(SSP-RK2/RK3) +- 局部阻力修正(入口/出口损失系数 ζ) +- HLLC、Roe、exact Riemann 等其他格式 +- 网格收敛性扫描 +- 真实气体状态方程 +- 多管网络、分叉、汇合 + +--- + +## 2. 物理与数值模型 + +### 2.1 0D 气瓶模型 + +对每个气瓶,假设内部状态均匀、气体静止,采用质量与总内能两个守恒量: + +``` +dm/dt = ṁ_in (质量守恒) +dU/dt = Ḣ_in (= ṁ·h_t,in) (能量守恒,开口系第一定律) +``` + +其中: +- `m` = 总质量 [kg] +- `U` = 总内能 [J],对理想气体 `U = P·V/(γ−1) = m·R·T/(γ−1)` +- `Ḣ` = 总焓流 [W],正值表示流入 + +派生量通过状态方程按需计算: + +``` +ρ = m/V +T = (U/m) · (γ−1)/R +P = ρ·R·T = U·(γ−1)/V +``` + +**设计取舍**:选 `(m, U)` 作底层状态而非 `(P, T)` 的理由见 §4.3。 + +### 2.2 1D 管路模型 + +一维可压缩欧拉方程(守恒形式): + +``` +∂W/∂t + ∂F(W)/∂x = 0 + +W = [ρ, ρu, ρE]ᵀ +F = [ρu, ρu² + P, u·(ρE + P)]ᵀ +``` + +其中 `E = e + u²/2`,`e = P/(ρ(γ−1))`。 + +离散化:有限体积法,cell-averaged piecewise constant(一阶重构)。 + +``` +W_i^{n+1} = W_i^n − (Δt/Δx)·(F_{i+1/2} − F_{i−1/2}) +``` + +- `N = 20` 个单元 +- `Δx = L/N = 0.05 m` +- 单元中心坐标 `x_i = (i + 0.5)·Δx`, i = 0..19 + +### 2.3 耦合:虚网格 + HLL + +**虚网格定义**:管路左右各增加一个"虚拟单元",每个时间步开头根据当前气瓶状态填充: + +``` +W_ghost_L = [ρ₁, 0, P₁/(γ−1)]ᵀ ← 高压瓶 +W_ghost_R = [ρ₂, 0, P₂/(γ−1)]ᵀ ← 低压瓶 +``` + +速度取 0 的假设:把气瓶视为滞止状态(stagnation state)的无穷大储气罐,真实的加速过程由 HLL 在界面上解出。这是**建模简化**,不是物理真相;其代价是瓶内到界面的速度-压力关系并非严格等熵,对短管、大压差场景误差可接受。 + +**注意区分两件事**: + +1. **系统总能守恒**:从"通量双用"的数学结构直接推出(见 §2.4),与 HLL 精度无关,严格到机器精度成立。 +2. **气瓶能量更新的物理一致性**:HLL 的 `flux[2] = u(ρE + P)` 等同于总焓流密度 `ρu·h_t`。对于 u_ghost=0 的滞止态气瓶,h_t,ghost = c_p·T_tank,这恰好是开口系能量方程对气瓶的正确源项——即便 Riemann 求解器引入数值耗散,能量守恒量也被 HLL 的守恒形式严格守住。 + +**HLL 数值通量**(Harten–Lax–van Leer): + +``` +S_L = min(u_L − a_L, u_R − a_R) +S_R = max(u_L + a_L, u_R + a_R) + + ┌ F_L if S_L ≥ 0 +F_HLL = │ (S_R·F_L − S_L·F_R + S_L·S_R·(W_R − W_L)) if S_L < 0 < S_R + │ ──────────────────────────────────── + │ S_R − S_L + └ F_R if S_R ≤ 0 +``` + +其中 `a = √(γP/ρ)`。 + +**关键性质**:HLL 自动处理双向流、激波、膨胀波和截流;无需手动判断流动方向或是否达到声速。 + +### 2.4 守恒性与通量双用 + +**核心耦合原理**:每个时间步,HLL 在管路左右两个边界各给出一个数值通量向量。**同一组通量被两边共享**: + +- 用于更新管路首/末单元的保守变量(有限体积更新) +- 乘以管路截面积 A 后,作为气瓶的 `ṁ, Ḣ` 源项更新气瓶 + +这一"通量双用"机制从数学上天然保证系统总质量和总能量守恒到机器精度,是验证实现正确性的唯一充分必要条件。 + +--- + +## 3. 时间推进算法 + +### 3.1 单步推进的 7 个阶段 + +``` +阶段 1 — 计算 CFL 时间步长 + a_i = √(γ P_i / ρ_i) for each cell i + a_max = max over i of (|u_i| + a_i) + Δt = CFL · Δx / a_max + Δt = min(Δt, t_end − t) 夹到 t_end + +阶段 2 — 冻结气瓶状态为虚网格 + W_ghost_L = tank1.ghost_state() # 本步用的高压瓶快照 + W_ghost_R = tank2.ghost_state() # 本步用的低压瓶快照 + (管路 W 的快照留给 pipe.step 内部处理,solver 不碰) + +阶段 3 — 计算两个边界通量(solver 层) + flux_L = HLL(W_ghost_L, pipe.W[:,0]) 左边界 (W_ghost_L ↔ cell 0) + flux_R = HLL(pipe.W[:,N-1], W_ghost_R) 右边界 (cell N-1 ↔ W_ghost_R) + + 注:只有这两个"边界通量"由 solver 计算,因为它们需要知道 ghost state。 + 所有 N-1 个"内部通量"由 pipe.step() 内部计算,solver 不参与。 + +阶段 4 — 管路一步推进(pipe.step 内部完成) + pipe.step(flux_L, flux_R, Δt): + W_snap = pipe.W.copy() + flux_int = [HLL(W_snap[:,i-1], W_snap[:,i]) for i in 1..N-1] 内部 N-1 个通量 + + 首单元: pipe.W[:,0] = W_snap[:,0] − (Δt/Δx)·(flux_int[1] − flux_L) + 末单元: pipe.W[:,N-1] = W_snap[:,N-1] − (Δt/Δx)·(flux_R − flux_int[N-1]) + 中间: pipe.W[:,i] = W_snap[:,i] − (Δt/Δx)·(flux_int[i+1] − flux_int[i]) + for i in 1..N-2 + +阶段 5 — 同一组边界通量更新气瓶 + fL_A = flux_L · A # (3,) 向量 × 标量面积 + fR_A = flux_R · A + tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=Δt, sign=-1) 流出 + tank2.apply_flux(mdot=fR_A[0], edot=fR_A[2], dt=Δt, sign=+1) 流入 + +阶段 6 — 推进时间 + t += Δt + +阶段 7 — 记录历史 + history['t'].append(t) + history['P1'].append(tank1.P) + history['T1'].append(tank1.T) + history['P2'].append(tank2.P) + history['T2'].append(tank2.T) + history['W_hist'].append(pipe.W.copy()) +``` + +### 3.2 关键顺序约束 + +1. **必须先一次性算完所有界面通量,再一次性更新所有单元**。如果边算边更新,会导致"新 W 和旧 W 混用"。 +2. **管路和气瓶的更新使用同一组冻结通量**。若用已被更新的 tank 状态重算左边界通量,即违反守恒性。 +3. **阶段 4 和阶段 5 的先后顺序不重要**(彼此独立),但不能交错。 + +### 3.3 CFL 估算与步数预估 + +初始时刻(管路全场 = 低压瓶态): +``` +a₀ = √(1.4 · 101325 / 1.177) ≈ 347 m/s +u₀ = 0 +Δt₀ = 0.5 · 0.05 / 347 ≈ 7.2 × 10⁻⁵ s +``` + +激波进入管路后 `|u| + a` 可达 600–800 m/s,Δt 收缩到 ~30–40 μs。 + +`t_end = 0.1 s` 预计对应 **约 2000–3000 步**。NumPy 向量化 N=20 的数组,预计全程 **<15 秒** 跑完。 + +--- + +## 4. 组件与模块划分 + +### 4.1 文件布局 + +``` +pipe-system-simulation-test/ +├── src/ +│ ├── config.py # 常量与工况(γ, R, V, P_init, L, D, N, CFL, t_end) +│ ├── riemann.py # HLL 数值通量(纯函数) +│ ├── tank.py # Tank 类(0D,状态 = (mass, U)) +│ ├── pipe.py # Pipe 类(1D 有限体积场 W, 推进方法) +│ ├── solver.py # 时间循环与气瓶-管路耦合编排 +│ ├── output.py # save_history / plot_timeseries / make_animation +│ └── main.py # 入口:装配、调度、sanity check +├── tests/ +│ ├── test_riemann.py +│ ├── test_tank.py +│ ├── test_pipe.py +│ └── test_integration.py +├── conftest.py # sys.path 注入,让 tests/ 能从 src/ 导入 +├── results/ # 运行产物(.npz, .png, .gif) +├── cases/ # (保留,后续扩展多工况用) +├── docs/ +│ └── superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md # 本文件 +└── scripts/ # (保留,后续工具脚本用) +``` + +### 4.2 模块职责 + +| 模块 | 职责 | 对外接口 | +|---|---|---| +| `config.py` | 纯数据:物理常数、工况、仿真控制参数 | 顶层常量 | +| `riemann.py` | HLL 数值通量计算 | `hll_flux(W_L, W_R, gamma) → ndarray(3,)` | +| `tank.py` | 0D 气瓶状态与演化 | `Tank` 类,`ghost_state()`, `apply_flux()`,派生属性 `P/T/rho` | +| `pipe.py` | 1D 管路状态与一步推进 | `Pipe` 类,`primitives()`, `max_wave_speed()`, `step()` | +| `solver.py` | 装配时间循环、调用 HLL、协调气瓶与管路更新 | `run(tank1, tank2, pipe, t_end, cfl, verbose=False, log_every=100) → history: dict` | +| `output.py` | 持久化与可视化 | `save_history()`, `plot_tank_timeseries()`, `make_pipe_animation()` | +| `main.py` | 入口脚本:构造对象 → 调 solver → 调 output → sanity check | `main()` | + +### 4.3 关键设计取舍 + +#### 4.3.1 Tank 底层状态选 `(mass, U)` 而非 `(P, T)` + +**理由**: + +- 守恒律在 `(m, U)` 空间是线性的:`m += ṁ·Δt`,`U += Ḣ·Δt`,两次加法即完成,零公式展开 +- `(P, T, ρ)` 三者由状态方程约束,只有 2 个自由度;若都当第一身份存,必须人工保证同步,任何一步漏更新就违反状态方程 +- HLL 的 `flux[2] = u(ρE + P)` 恰好等于总焓流密度 `ρu·h_t`,乘以 A 后就是开口系能量方程的右端项 `Ḣ`;写成 `dU/dt = Ḣ` 直接对应代码 `tank.U += flux[2]·A·dt`,无需手工添加"流动功修正" +- 派生属性(`tank.P`, `tank.T`, `tank.rho`)通过 `@property` 实时计算,永远与 `(m, U)` 一致,外部读不到过时值 + +**代价**:每次访问 `tank.P` 需要一次除法 + 一次乘法。对 0D 气瓶完全可忽略。 + +#### 4.3.2 界面通量"双用",由 solver 显式编排 + +- `solver.run()` 在每步先算边界 flux,再把 flux 分发给 `pipe.step(flux_L, flux_R, dt)` 和 `tank.apply_flux(ṁ, Ḣ, dt, sign)` +- `pipe.step()` 内部负责所有**内部**界面通量的计算和管路单元更新 +- 气瓶不知道管路存在,管路不知道气瓶存在;耦合知识集中在 solver 一处 + +#### 4.3.3 ghost_state() 作为 Tank 的方法 + +- "如何把自身暴露成虚网格"是 tank 的职责,不是 solver 的职责 +- solver 代码里就是 `hll_flux(tank1.ghost_state(), pipe.W[:,0], γ)`,可读性最高 +- 未来如果要扩展成"带局部阻力修正的虚网格",只改 `Tank.ghost_state()` 即可 + +#### 4.3.4 sign 参数约定 + +`Tank.apply_flux(mdot, edot, dt, sign)` 接受一个显式的 sign 参数: + +- `sign = -1`:气流"流出"瓶子。对 tank1(高压瓶),左边界通量向右为正,即流出,所以用 -1 +- `sign = +1`:气流"流入"瓶子。对 tank2(低压瓶),右边界通量向右为正,即流入,所以用 +1 + +显式 sign 比在 Tank 内部判断"我是 tank1 还是 tank2"更干净。 + +--- + +## 5. 数据流与历史记录 + +### 5.1 History 数据结构 + +```python +history = { + 't': np.ndarray, # shape (n_steps,) 时间序列 + 'P1': np.ndarray, # shape (n_steps,) 高压瓶压力 + 'T1': np.ndarray, # shape (n_steps,) 高压瓶温度 + 'P2': np.ndarray, # shape (n_steps,) 低压瓶压力 + 'T2': np.ndarray, # shape (n_steps,) 低压瓶温度 + 'W_hist': np.ndarray, # shape (n_steps, 3, N) 管路全程守恒变量 +} +``` + +循环中用 Python list + append 收集,循环结束后一次性 `np.stack` / `np.asarray` 转成数组。 + +**内存估算**: +``` +n_steps ≈ 3000, 3 vars × 20 cells × 8 bytes +→ 每步 480 B, 全程 ≈ 1.4 MB +``` +完全无压力,**每步都记**,不做降采样。 + +### 5.2 持久化:`results/history.npz` + +```python +np.savez_compressed( + "results/history.npz", + t = history['t'], + P1 = history['P1'], T1 = history['T1'], + P2 = history['P2'], T2 = history['T2'], + W_hist = history['W_hist'], + x = pipe_cell_centers, # shape (N,) + dx = pipe.dx, + gamma = GAMMA, + R_gas = R_GAS, + area = pipe.area, +) +``` + +未来加载与访问示例: +```python +d = np.load("results/history.npz") +rho = d['W_hist'][:, 0, :] # (n_steps, N) +u = d['W_hist'][:, 1, :] / rho +P = (d['W_hist'][:, 2, :] - 0.5*rho*u**2) * (d['gamma'] - 1) +``` + +### 5.3 可视化产物 + +``` +results/ +├── history.npz 全量时间序列,可离线复盘 +├── tank_pressure.png P1(t), P2(t) 双曲线 +├── tank_temperature.png T1(t), T2(t) 双曲线 +└── pipe_animation.gif 管内 P(x), u(x), T(x) 三子图演化 +``` + +**动画细节**: + +- 用 `matplotlib.animation.FuncAnimation` + `PillowWriter`(GIF 输出,免 ffmpeg 依赖) +- 布局:3 个纵向子图,共享 x 轴(管路空间坐标 0–1 m) +- 每帧标题:`t = X.XXXX s (step k/n_steps)` +- 默认 `stride = 10`,约 300 帧,GIF 文件预计 5–10 MB +- 数据存储仍为全量;stride 只影响动画帧数 + +--- + +## 6. 错误处理 + +### 6.1 验证点清单 + +| 位置 | 检查 | 失败时 | +|---|---|---| +| 程序启动 | `results/` 目录存在(`os.makedirs(exist_ok=True)`) | 静默创建 | +| `config.py` 顶层 | γ>1, R>0, V>0, L>0, D>0, N≥2, P>0, T>0, CFL∈(0,1] | `AssertionError` | +| 时间步开始 | 管路所有单元 ρ>0, P>0 | `RuntimeError("负密度/负压力 at step k")` | +| CFL 步长 | `Δt > 1e-12` | `RuntimeError("dt 退化")` | +| HLL 内部 | 恢复原始变量时 ρ>0, P>0 | `ValueError` + 打印 W_L, W_R | +| Tank 更新后 | `mass > 0` | `RuntimeError("气瓶质量非正")` | +| 仿真结束 | `n_steps > 0` | `RuntimeError("一步都没跑")` | +| main.py 末尾 | 总质量相对误差 < 1e-10 | `AssertionError`(作为 sanity check 打印) | +| main.py 末尾 | 总能量相对误差 < 1e-10 | `AssertionError`(作为 sanity check 打印) | + +### 6.2 策略 + +- **异常直接往上抛**,不 catch、不重试、不降级。让 Python traceback 定位问题 +- **不** `try/except` 包装 `main.py` +- `solver.run()` 提供 `verbose=False` 开关(默认关)和 `log_every=100` 打印间隔,用于调试时查看步进日志 + +### 6.3 明确不做的检查 + +- 负质量流量判断(HLL 自动处理双向) +- 显式判截流(Riemann 自动处理) +- 动态 CFL 调整策略 +- 并行/线程安全 + +--- + +## 7. 测试策略 + +### 7.1 测试清单(共 10 条,4 个文件,预计 <5 秒) + +#### `tests/test_riemann.py` — HLL 求解器 + +1. **左右状态相同 → 返回纯物理通量**(无数值耗散) +2. **静止接触间断**(两侧 u=0,仅 ρ 不同)→ 质量通量和能量通量应为 0,动量通量 = P +3. **Sod 激波管初值的符号与量级检查**(不做精确值断言,检查 flux[0], flux[1], flux[2] 均 > 0) + +#### `tests/test_tank.py` — 气瓶 + +4. **初始状态自洽**:给定 (P, T, V) 构造 → `mass` 和 `U` 满足理想气体关系 +5. **ghost_state 是滞止态**:动量项 = 0,能量项 = P/(γ−1) +6. **apply_flux 方向性**:sign=-1 应减质量、减内能、降压 + +#### `tests/test_pipe.py` — 管路 + +7. **均匀初始化**:所有单元的 W 完全相同;`primitives()` 返回 u=0, P=P_init +8. **静止状态 + 一致压力通量 → 一步后 W 不变**(机器精度) + +#### `tests/test_integration.py` — 端到端守恒律 + +9. **短仿真总质量守恒**:`t_end=1e-3`,`|Δm_total/m_total| < 1e-10` +10. **短仿真总能量守恒**:同上,`|ΔU_total/U_total| < 1e-10` + +### 7.2 守恒律断言的价值 + +守恒律是**数学真理**,任何正确实现都必然满足。如果实现中"通量双用"机制写错(例如 tank 用了被更新后的 flux,或算 flux 时没快照 W),守恒律立刻以可观察的数量级被破坏。这比"和参考解对比"更可靠——后者有过拟合风险。 + +### 7.3 不做的测试 + +- 绘图/动画输出(肉眼验证更高效) +- `main.py` 装配逻辑(无独立逻辑) +- 网格收敛 / CFL 收敛 / 阶数验证 +- Method of Manufactured Solutions +- mock/fixture 框架(10 条测试不需要抽象) + +### 7.4 运行方式 + +```bash +pytest tests/ -v +``` + +项目根运行。`conftest.py` 自动注入 `src/` 到 `sys.path`。 + +### 7.5 可测性要求(对生产代码的约束) + +1. 各模块顶层**无副作用**(禁止 print/read-file/random-seed) +2. `solver.run()` 不读取 `config` 模块,所有参数通过入参传入 +3. `Pipe.primitives()` 返回 `(ρ, u, P, a)` 四元组,方便测试断言 + +--- + +## 8. 验收标准 + +MVP 完成的判据: + +1. ✅ `python src/main.py` 能跑完、不崩溃、在 `results/` 下生成全部 4 个文件 +2. ✅ `pytest tests/ -v` 全部 10 条测试通过 +3. ✅ main.py 末尾的"总质量 & 总能量相对误差 < 1e-10" sanity check 通过 +4. ✅ `P1(t)` 单调递减、`P2(t)` 单调递增 +5. ✅ 动画能用常见 GIF 播放器(浏览器、图片查看器)正常打开 +6. ✅ `history.npz` 可以用 `np.load` 重新加载并访问所有字段 +7. ✅ 全程运行时间(不含测试) < 30 秒 + +--- + +## 9. 已知局限与后续扩展 + +**本 MVP 的物理局限**: + +- 细管(D=5mm)下摩擦损失可能显著,当前绝热无摩擦会高估压力传递速率 +- 瓶内气体从滞止到出口的加速不满足严格等熵(因为 u_ghost=0 的简化) +- 一阶空间导致膨胀波被耗散 +- 一阶时间精度对大 CFL 下的波传播有额外耗散 + +**后续可扩展方向**(不在本 MVP 范围内): + +1. 加 Darcy–Weisbach 摩擦源项(`src/friction.py`) +2. 加壁面对流换热(`src/heat_transfer.py`) +3. 升级到 MUSCL + minmod 二阶空间 +4. 升级到 SSP-RK2 二阶时间 +5. 在 `Tank.ghost_state()` 中加入入口损失系数 ζ +6. 配置化工况:YAML/JSON 读入 `cases/*.yaml` +7. 扩展到多管网络(带分叉节点的 0D/1D 混合拓扑) +8. 真实气体状态方程(Peng–Robinson / GERG-2008) + +--- + +## 10. 参考 + +- Toro, E. F. *Riemann Solvers and Numerical Methods for Fluid Dynamics*. Springer. (HLL 格式的标准参考) +- Harten, A., Lax, P. D., and van Leer, B. (1983). "On upstream differencing and Godunov-type schemes for hyperbolic conservation laws." *SIAM Review* 25(1): 35–61. +- LeVeque, R. J. *Finite Volume Methods for Hyperbolic Problems*. Cambridge University Press. +- 用户在 brainstorming 阶段提供的参考资料(涵盖 0D–1D 耦合的虚网格 + Riemann 方法) diff --git a/docs/superpowers/specs/2026-04-16-cryo-tank-module-design.md b/docs/superpowers/specs/2026-04-16-cryo-tank-module-design.md new file mode 100644 index 0000000..8c1bf3e --- /dev/null +++ b/docs/superpowers/specs/2026-04-16-cryo-tank-module-design.md @@ -0,0 +1,401 @@ +# 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 diff --git a/docs/管路系统仿真技术说明文档.docx b/docs/管路系统仿真技术说明文档.docx new file mode 100644 index 0000000..d72b180 Binary files /dev/null and b/docs/管路系统仿真技术说明文档.docx differ diff --git a/generate_doc.py b/generate_doc.py new file mode 100644 index 0000000..a6a1fb4 --- /dev/null +++ b/generate_doc.py @@ -0,0 +1,821 @@ +#!/usr/bin/env python3 +"""Generate a Word document describing the pipe system simulation theory and implementation.""" + +import os +import sys + +from docx import Document +from docx.shared import Pt, Inches, Cm, RGBColor +from docx.enum.text import WD_ALIGN_PARAGRAPH +from docx.enum.table import WD_TABLE_ALIGNMENT +from docx.oxml.ns import qn + + +def set_cell_shading(cell, color_hex): + shading_elm = cell._element.get_or_add_tcPr() + shd = shading_elm.makeelement(qn('w:shd'), { + qn('w:val'): 'clear', + qn('w:color'): 'auto', + qn('w:fill'): color_hex, + }) + shading_elm.append(shd) + + +def add_equation(doc, text, bold=False): + """Add a centered equation paragraph with Cambria Math font.""" + p = doc.add_paragraph() + p.alignment = WD_ALIGN_PARAGRAPH.CENTER + p.paragraph_format.space_before = Pt(6) + p.paragraph_format.space_after = Pt(6) + run = p.add_run(text) + run.font.name = 'Cambria Math' + run.font.size = Pt(11) + if bold: + run.bold = True + return p + + +def add_code_block(doc, code): + """Add a code block with monospace font and gray background.""" + p = doc.add_paragraph() + p.paragraph_format.space_before = Pt(4) + p.paragraph_format.space_after = Pt(4) + p.paragraph_format.left_indent = Cm(1) + run = p.add_run(code) + run.font.name = 'Consolas' + run.font.size = Pt(9) + run.font.color.rgb = RGBColor(0x33, 0x33, 0x33) + return p + + +def make_table(doc, headers, rows, col_widths=None): + """Create a formatted table.""" + table = doc.add_table(rows=1 + len(rows), cols=len(headers)) + table.style = 'Table Grid' + table.alignment = WD_TABLE_ALIGNMENT.CENTER + + # Header row + for j, h in enumerate(headers): + cell = table.rows[0].cells[j] + cell.text = h + set_cell_shading(cell, 'D9E2F3') + for paragraph in cell.paragraphs: + for run in paragraph.runs: + run.bold = True + run.font.size = Pt(10) + + # Data rows + for i, row in enumerate(rows): + for j, val in enumerate(row): + cell = table.rows[i + 1].cells[j] + cell.text = str(val) + for paragraph in cell.paragraphs: + for run in paragraph.runs: + run.font.size = Pt(10) + + if col_widths: + for i, w in enumerate(col_widths): + for row in table.rows: + row.cells[i].width = Cm(w) + + return table + + +def build_document(): + doc = Document() + + # --- Global style --- + style = doc.styles['Normal'] + style.font.name = 'Times New Roman' + style.font.size = Pt(11) + style.paragraph_format.line_spacing = 1.15 + style.paragraph_format.space_after = Pt(6) + + for level in range(1, 4): + hs = doc.styles[f'Heading {level}'] + hs.font.name = 'Times New Roman' + hs.font.color.rgb = RGBColor(0x1F, 0x3A, 0x5F) + + # ======================================================================== + # TITLE + # ======================================================================== + title = doc.add_heading('0D-1D Tank-Pipe Coupled Transient Simulation\nTechnical Documentation', level=0) + title.alignment = WD_ALIGN_PARAGRAPH.CENTER + for run in title.runs: + run.font.size = Pt(22) + + p = doc.add_paragraph() + p.alignment = WD_ALIGN_PARAGRAPH.CENTER + run = p.add_run('0D-1D Tank-Pipe Blowdown Simulation\nTheory, Numerical Methods, and Source Code Description') + run.font.size = Pt(12) + run.font.color.rgb = RGBColor(0x66, 0x66, 0x66) + run.italic = True + + doc.add_page_break() + + # ======================================================================== + # 1. SYSTEM OVERVIEW + # ======================================================================== + doc.add_heading('1 System Overview', level=1) + + doc.add_paragraph( + 'This document describes the theoretical foundations, numerical schemes, ' + 'and source code implementation of a 0D-1D coupled transient simulation system. ' + 'The system models the blowdown process where gas flows from a high-pressure ' + 'tank (Tank 1) through a slender pipe into a low-pressure tank (Tank 2).' + ) + + doc.add_heading('1.1 Physical Scenario', level=2) + doc.add_paragraph( + 'The system consists of three components:\n' + '(1) High-pressure tank (upstream, Tank 1): volume V1 = 5 m^3, initial pressure P1 = 5 MPa, temperature T1 = 300 K.\n' + '(2) Pipe: length L = 1 m, inner diameter D = 5 mm, divided into N = 20 finite-volume cells.\n' + '(3) Low-pressure tank (downstream, Tank 2): volume V2 = 10 m^3, initial pressure P2 = 2 MPa, temperature T2 = 300 K.\n' + '\n' + 'The two tanks are connected through the pipe. Due to the pressure difference, ' + 'gas flows from the high-pressure side to the low-pressure side. ' + 'The flow inside the pipe is one-dimensional compressible flow, while the two tanks, ' + 'whose volumes are much larger than the pipe, have spatially uniform internal states ' + 'and are modeled with 0D lumped-parameter models.' + ) + + doc.add_heading('1.2 Modeling Strategy', level=2) + doc.add_paragraph( + 'The overall approach is "0D-1D coupling":\n' + '- Tanks: 0D Lumped Parameter Model -- internal state (P, T, rho) is spatially uniform, ' + 'varies only with time.\n' + '- Pipe: 1D Compressible Euler Equations -- captures pressure waves, shock waves, ' + 'and expansion waves propagating through the pipe.\n' + '- Coupling: Ghost Cell Method -- a virtual cell is placed at each end of the pipe, ' + 'filled with the current tank state, then the HLL Riemann solver computes boundary fluxes.\n' + '- Friction: Source Term Method (Operator Splitting) -- Darcy-Weisbach wall friction ' + 'added to the momentum equation.' + ) + + doc.add_heading('1.3 Code Module Structure', level=2) + + make_table(doc, + ['Module', 'Responsibility', 'Core Class/Function'], + [ + ['config.py', 'Physical constants, geometry, simulation control', '14 constants + assertion checks'], + ['tank.py', '0D tank model', 'Tank class (mass, U, ghost_state, apply_flux)'], + ['riemann.py', 'HLL Riemann numerical flux', 'hll_flux(W_L, W_R, gamma)'], + ['pipe.py', '1D pipe finite-volume model', 'Pipe class (W, primitives, step)'], + ['friction.py', 'Darcy-Weisbach friction factor', 'darcy_friction_factor(Re, eps_D)'], + ['solver.py', 'Time-stepping driver (coupling orchestrator)', 'run(tank1, tank2, pipe, ...)'], + ['output.py', 'Post-processing: storage, plots, animation, report', 'save_history, plot_*, make_pipe_animation'], + ['main.py', 'Entry point: assemble, run, verify, output', 'main()'], + ], + col_widths=[3.5, 6, 6.5], + ) + + doc.add_page_break() + + # ======================================================================== + # 2. 0D TANK MODEL + # ======================================================================== + doc.add_heading('2 0D Tank Model', level=1) + + doc.add_heading('2.1 Basic Assumptions', level=2) + doc.add_paragraph( + 'The tanks use a lumped-parameter (0D) model with the following assumptions:\n' + '(1) The gas is an ideal gas with equation of state P = rho * R * T.\n' + '(2) Internal state is spatially uniform -- density rho, temperature T, and pressure P ' + 'are functions of time only.\n' + '(3) The macroscopic velocity inside the tank is zero (u = 0); kinetic energy is negligible.\n' + '(4) The tank walls are adiabatic; no heat conduction to the environment.\n' + '(5) The tank volume is fixed.' + ) + + doc.add_heading('2.2 Conservation Equations', level=2) + doc.add_paragraph( + 'For an open-system tank, mass conservation and energy conservation ' + '(open-system first law of thermodynamics) are:' + ) + add_equation(doc, 'dm/dt = m_dot_in') + add_equation(doc, 'dU/dt = H_dot_in = m_dot_in * h_t,in') + doc.add_paragraph( + 'Where:\n' + '- m = rho * V is the total gas mass in the tank [kg]\n' + '- U = P*V/(gamma-1) is the total internal energy [J] (since u=0, total energy = internal energy)\n' + '- m_dot_in is the mass flow rate entering the tank [kg/s]\n' + '- H_dot_in is the total enthalpy flow rate entering the tank [W]\n' + '- h_t,in is the specific total enthalpy of the incoming gas [J/kg]\n' + '\n' + 'Key design: In the code, the Tank primary state is (mass, U), ' + 'while rho, T, P are all computed on demand via @property:' + ) + add_equation(doc, 'rho = mass / V') + add_equation(doc, 'T = (U / mass) * (gamma - 1) / R') + add_equation(doc, 'P = rho * R * T') + doc.add_paragraph( + 'This design avoids state-synchronization bugs after time-step updates -- ' + 'only mass and U are updated; all derived quantities are automatically consistent.' + ) + + doc.add_heading('2.3 Code Implementation (tank.py)', level=2) + doc.add_paragraph('The Tank constructor computes initial mass and U from (V, P_init, T_init, gamma, R):') + add_code_block(doc, + 'rho = P_init / (R_gas * T_init)\n' + 'self.mass = rho * V\n' + 'self.U = P_init * V / (gamma - 1)') + + doc.add_paragraph( + 'The apply_flux(mdot, edot, dt, sign) method performs each time-step update:\n' + ' mass += sign * mdot * dt\n' + ' U += sign * edot * dt\n' + 'Where sign = -1 means gas flows out (upstream Tank 1), ' + 'sign = +1 means gas flows in (downstream Tank 2).' + ) + + doc.add_page_break() + + # ======================================================================== + # 3. 1D PIPE MODEL + # ======================================================================== + doc.add_heading('3 1D Pipe Model', level=1) + + doc.add_heading('3.1 Governing Equations: 1D Compressible Euler Equations', level=2) + doc.add_paragraph( + 'Gas flow inside the pipe is described by the 1D compressible Euler equations ' + '(without friction):' + ) + add_equation(doc, 'dW/dt + dF(W)/dx = 0') + doc.add_paragraph('Where the conservative variable vector W and flux vector F(W) are:') + add_equation(doc, 'W = [rho, rho*u, rho*E]^T') + add_equation(doc, 'F = [rho*u, rho*u^2 + P, u*(rho*E + P)]^T') + doc.add_paragraph( + 'Physical meaning of each component:\n' + '- W[0] = rho: density [kg/m^3]\n' + '- W[1] = rho*u: momentum density [kg/(m^2*s)]\n' + '- W[2] = rho*E: total energy density [J/m^3], where E = e + u^2/2 is specific total energy, ' + 'e = P/(rho*(gamma-1)) is specific internal energy\n' + '\n' + 'The three components of F correspond to mass flux, momentum flux (including pressure), ' + 'and energy flux, respectively.' + ) + doc.add_paragraph( + 'Equation of state (ideal gas closure):' + ) + add_equation(doc, 'P = (gamma - 1) * (rho*E - 0.5*rho*u^2)') + add_equation(doc, 'a = sqrt(gamma * P / rho) [speed of sound]') + + doc.add_heading('3.2 Finite Volume Discretization', level=2) + doc.add_paragraph( + 'The pipe is uniformly divided into N finite-volume cells, each of width dx = L/N. ' + 'Within each cell, the conservative variables W take cell-averaged values ' + '(piecewise-constant reconstruction), giving first-order spatial accuracy.' + ) + doc.add_paragraph( + 'For the i-th cell [x_{i-1/2}, x_{i+1/2}], integrating the conservation law ' + 'yields the semi-discrete form:' + ) + add_equation(doc, 'dW_i/dt = -(1/dx) * [F_{i+1/2} - F_{i-1/2}]') + doc.add_paragraph( + 'Where F_{i+1/2} is the numerical flux at the interface between cells i and i+1. ' + 'The key challenge is that adjacent cell states are generally discontinuous at the interface ' + '(a Riemann discontinuity), and a Riemann solver is needed to determine the interface flux.' + ) + + doc.add_heading('3.3 Code Implementation (pipe.py)', level=2) + doc.add_paragraph( + 'The Pipe class stores the W array with shape (3, N), where each column is one cell\'s ' + 'conservative variables. The initial state is computed from (P_init, T_init):' + ) + add_code_block(doc, + 'rho = P_init / (R_gas * T_init)\n' + 'E_density = P_init / (gamma - 1) # u=0, total energy density = internal\n' + 'W[0, :] = rho # density\n' + 'W[1, :] = 0.0 # momentum density (u=0)\n' + 'W[2, :] = E_density # total energy density') + doc.add_paragraph( + 'The primitives() method recovers primitive variables (rho, u, P, a) from W. ' + 'The step(flux_L, flux_R, dt) method performs one time-step advancement.' + ) + + doc.add_page_break() + + # ======================================================================== + # 4. HLL RIEMANN SOLVER + # ======================================================================== + doc.add_heading('4 HLL Riemann Solver', level=1) + + doc.add_heading('4.1 The Riemann Problem', level=2) + doc.add_paragraph( + 'In the finite-volume method, the states on the left and right sides of each cell interface ' + 'are generally different, forming a local Riemann problem. For the 1D Euler equations, ' + 'the exact Riemann solution contains three waves: a left-going shock/rarefaction, ' + 'a contact discontinuity, and a right-going shock/rarefaction.' + '\n\n' + 'Exact Riemann solvers are computationally expensive, so engineering practice typically uses ' + 'approximate Riemann solvers. This system uses the HLL (Harten-Lax-van Leer) scheme, ' + 'which is simple and robust.' + ) + + doc.add_heading('4.2 HLL Scheme Derivation', level=2) + doc.add_paragraph( + 'The HLL scheme assumes the Riemann fan is separated by two waves (S_L and S_R) ' + 'into three regions:\n' + '- x/t < S_L: left state W_L (undisturbed region ahead of waves)\n' + '- S_L < x/t < S_R: intermediate state W* ("HLL average state")\n' + '- x/t > S_R: right state W_R (undisturbed region ahead of waves)\n' + '\n' + 'The flux at the interface x=0 depends on the signs of S_L and S_R:' + ) + + doc.add_paragraph('Case 1: If S_L >= 0 (both waves travel rightward), the interface is left of all waves:') + add_equation(doc, 'F_HLL = F(W_L)') + + doc.add_paragraph('Case 2: If S_R <= 0 (both waves travel leftward), the interface is right of all waves:') + add_equation(doc, 'F_HLL = F(W_R)') + + doc.add_paragraph('Case 3: If S_L < 0 < S_R (interface is between the two waves), weighted average:') + add_equation(doc, 'F_HLL = [S_R*F_L - S_L*F_R + S_L*S_R*(W_R - W_L)] / (S_R - S_L)') + + doc.add_heading('4.3 Wave Speed Estimates', level=2) + doc.add_paragraph( + 'The choice of wave speeds S_L and S_R is critical to the HLL scheme. ' + 'This system uses the Davis estimate:' + ) + add_equation(doc, 'S_L = min(u_L - a_L, u_R - a_R)') + add_equation(doc, 'S_R = max(u_L + a_L, u_R + a_R)') + doc.add_paragraph( + 'Where u is the flow velocity and a = sqrt(gamma*P/rho) is the local speed of sound. ' + 'This estimate is simple and robust, ensuring all physical signal propagation speeds ' + 'are contained within [S_L, S_R].' + ) + + doc.add_heading('4.4 Code Implementation (riemann.py)', level=2) + doc.add_paragraph('The hll_flux(W_L, W_R, gamma) function follows this procedure:') + doc.add_paragraph( + '(1) Recover primitive variables (rho, u, P, a) from W_L and W_R\n' + '(2) Compute physical fluxes F_L = F(W_L), F_R = F(W_R)\n' + '(3) Estimate wave speeds S_L, S_R\n' + '(4) Return the HLL flux according to the three cases above\n' + '(5) Input validation: raise ValueError if rho <= 0 or P <= 0' + ) + + doc.add_page_break() + + # ======================================================================== + # 5. GHOST CELL BOUNDARY -- THE KEY COUPLING MECHANISM + # ======================================================================== + doc.add_heading('5 Boundary Conditions: Ghost Cell Method', level=1) + + doc.add_paragraph( + 'This is one of the most critical design elements of the simulation: how to couple ' + 'the 0D tank models with the 1D pipe model. The approach used here is the ' + 'Ghost Cell Method. The basic idea is to place a "virtual cell" at each end of the pipe, ' + 'fill it with the current tank state, and then call the Riemann solver between the ' + 'virtual cell and the first/last real pipe cell to compute the boundary numerical flux.' + ) + + doc.add_heading('5.1 Ghost State Construction', level=2) + doc.add_paragraph( + 'For a tank connected to the pipe, the Ghost State is generated by ' + 'Tank.ghost_state(). The conservative variables of the ghost cell are:' + ) + add_equation(doc, 'W_ghost = [rho_tank, 0, P_tank / (gamma-1)]^T') + doc.add_paragraph( + 'Where:\n' + '- rho_tank is the current tank density\n' + '- Momentum rho*u = 0 (stagnation assumption: macroscopic velocity inside the tank is zero)\n' + '- rho*E = P/(gamma-1) (since u=0, total energy density = internal energy density)\n' + '\n' + 'The "Stagnation Assumption" is the key here:\n' + 'The tank volume is much larger than the pipe cross-section times pipe length, ' + 'so the bulk gas velocity inside the tank is extremely low -- the tank side of the pipe ' + 'inlet can be treated as a stagnation state. ' + 'The Riemann solver automatically determines the correct flow direction and flux magnitude ' + 'based on the state difference between the pipe side and the tank side.' + ) + + doc.add_heading('5.2 Left Boundary (Tank 1 -> Pipe Inlet)', level=2) + doc.add_paragraph('The specific computation steps:') + doc.add_paragraph( + '(1) Get ghost state from Tank 1: W_ghost_L = tank1.ghost_state()\n' + '(2) First real pipe cell state: W_pipe_0 = pipe.W[:, 0]\n' + '(3) Call HLL solver: flux_L = hll_flux(W_ghost_L, W_pipe_0, gamma)\n' + '(4) Here HLL "left state" = tank ghost cell, "right state" = pipe cell 0\n' + '\n' + 'The returned flux_L is a 3-component vector [mass flux, momentum flux, energy flux], ' + 'with positive direction from left to right (i.e., from tank into pipe).' + ) + + p = doc.add_paragraph() + run = p.add_run( + 'Physical interpretation: Since Tank 1 pressure (5 MPa) is much higher than ' + 'the pipe initial pressure (2 MPa), the Riemann solver automatically produces ' + 'a positive left-to-right flux, driving gas from the high-pressure tank into the pipe.' + ) + run.italic = True + + doc.add_heading('5.3 Right Boundary (Pipe Outlet -> Tank 2)', level=2) + doc.add_paragraph( + '(1) Last real pipe cell state: W_pipe_{N-1} = pipe.W[:, -1]\n' + '(2) Get ghost state from Tank 2: W_ghost_R = tank2.ghost_state()\n' + '(3) Call HLL solver: flux_R = hll_flux(W_pipe_{N-1}, W_ghost_R, gamma)\n' + '(4) Here HLL "left state" = last pipe cell, "right state" = tank ghost cell\n' + '\n' + 'flux_R positive direction is also left to right (i.e., from pipe into tank).' + ) + + doc.add_heading('5.4 Schematic Diagram', level=2) + doc.add_paragraph( + 'The following diagram shows the spatial layout of the ghost cell method:' + ) + add_code_block(doc, + ' [Tank 1] | Cell 0 | Cell 1 | ... | Cell N-1 | [Tank 2]\n' + ' (ghost_L) | | | | | (ghost_R)\n' + ' ^ ^\n' + ' flux_L flux_R\n' + ' = hll(ghost_L, = hll(pipe[:,-1],\n' + ' pipe[:,0]) ghost_R)\n' + ) + doc.add_paragraph( + 'Note: The ghost cells do not occupy physical space. They only provide ' + '"outside-the-pipe" information for the HLL solver to correctly compute boundary fluxes. ' + 'The ghost cell W values are re-frozen (snapshot) at the start of each time step ' + 'to ensure boundary conditions remain consistent within a single time step.' + ) + + doc.add_page_break() + + # ======================================================================== + # 6. FLUX DOUBLING + # ======================================================================== + doc.add_heading('6 Flux Doubling: Guaranteeing Conservation', level=1) + + doc.add_paragraph( + 'This is the core mechanism ensuring physical conservation laws ' + '(mass conservation, energy conservation) in the entire simulation system.' + ) + + doc.add_heading('6.1 The Problem', level=2) + doc.add_paragraph( + 'In 0D-1D coupling, the pipe boundary flux must be used not only to update ' + 'the pipe internal cell states, but also to simultaneously update the adjacent ' + 'tank\'s mass and energy. If the pipe and the tank use different flux calculations ' + 'to update their own states, then the total system mass and energy cannot be exactly ' + 'conserved -- an unphysical phenomenon of "mass appearing or disappearing" ' + 'would occur at the interface.' + ) + + doc.add_heading('6.2 Solution: Same Flux Updates Both Sides', level=2) + doc.add_paragraph( + 'The approach in this system is: compute the boundary HLL flux only once, ' + 'then use the same flux vector to simultaneously update:\n' + '(1) The pipe boundary cell via the finite-volume formula (in pipe.step using flux_L, flux_R)\n' + '(2) The tank mass and energy (in tank.apply_flux)\n' + '\n' + 'This is called "Flux Doubling".' + ) + + doc.add_paragraph('Using the left boundary as an example, the key logic in solver.py:') + add_code_block(doc, + '# Phase 3: compute left boundary flux (only once)\n' + 'flux_L = hll_flux(W_ghost_L, pipe.W[:, 0], gamma)\n' + '\n' + '# Phase 4: pipe uses this flux for update (as left interface flux)\n' + 'pipe.step(flux_L, flux_R, dt)\n' + ' -> inside: W[:,0] -= (dt/dx) * (flux_int[:,0] - flux_L)\n' + '\n' + '# Phase 5: tank also uses the SAME flux for update\n' + 'fL_A = flux_L * pipe.area # flux x cross-section area = physical rate\n' + 'tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=dt, sign=-1)\n' + ' -> inside: mass -= fL_A[0] * dt\n' + ' U -= fL_A[2] * dt') + + doc.add_heading('6.3 Why Does This Guarantee Conservation?', level=2) + doc.add_paragraph( + 'Consider mass conservation at the left boundary. In one time step dt:\n' + '- Pipe cell 0 gains mass from left interface flux = +flux_L[0] * A * dt\n' + '- Tank 1 loses mass from the same flux = flux_L[0] * A * dt\n' + '\n' + 'The two cancel exactly. The same applies to the right boundary and the energy component. ' + 'Therefore, the total system quantity (Tank 1 + all Pipe cells + Tank 2) is strictly ' + 'conserved at every time step, with errors only from floating-point roundoff, ' + 'typically at the 10^(-15) level.' + ) + + doc.add_paragraph( + 'Verification result: After simulation, the relative errors of total mass and total ' + 'energy are ~6.5 x 10^(-16) (mass) and ~1.1 x 10^(-15) (energy), i.e., ' + 'machine-epsilon level.' + ) + + doc.add_page_break() + + # ======================================================================== + # 7. TIME STEPPING + # ======================================================================== + doc.add_heading('7 Time-Stepping Algorithm', level=1) + + doc.add_heading('7.1 Explicit Euler Time Integration', level=2) + doc.add_paragraph( + 'This system uses first-order explicit Euler time integration. For each pipe cell:' + ) + add_equation(doc, 'W_i^(n+1) = W_i^(n) - (dt/dx) * [F_{i+1/2}^(n) - F_{i-1/2}^(n)]') + doc.add_paragraph( + 'For the tanks:' + ) + add_equation(doc, 'mass^(n+1) = mass^(n) + sign * m_dot * dt') + add_equation(doc, 'U^(n+1) = U^(n) + sign * E_dot * dt') + + doc.add_heading('7.2 CFL Condition', level=2) + doc.add_paragraph( + 'Explicit time integration has a stability restriction: the time step cannot exceed ' + 'the time for a signal to traverse one grid cell. ' + 'The CFL (Courant-Friedrichs-Lewy) condition is:' + ) + add_equation(doc, 'dt = CFL * dx / max_i(|u_i| + a_i)') + doc.add_paragraph( + 'Where CFL is a safety factor in (0, 1] (this system uses CFL = 0.5), ' + 'and max_i(|u_i| + a_i) is the maximum characteristic speed across all pipe cells. ' + 'Additionally, dt is clamped to ensure it does not exceed the remaining simulation time ' + 't_end - t.' + ) + + doc.add_heading('7.3 Complete 7-Phase Time Step', level=2) + doc.add_paragraph( + 'Each time step follows this execution flow (corresponding to the run() function in solver.py):' + ) + + make_table(doc, + ['Phase', 'Operation', 'Description'], + [ + ['Phase 1', 'CFL dt calculation', 'dt = CFL*dx/max(|u|+a), clamped to t_end-t'], + ['Phase 2', 'Freeze ghost cells', 'W_ghost_L = tank1.ghost_state()\nW_ghost_R = tank2.ghost_state()'], + ['Phase 3', 'Compute boundary HLL fluxes', 'flux_L = hll(ghost_L, pipe[:,0])\nflux_R = hll(pipe[:,-1], ghost_R)'], + ['Phase 4', 'Advance pipe', 'pipe.step(flux_L, flux_R, dt)\nincludes internal HLL fluxes + friction source'], + ['Phase 5', 'Advance tanks (flux doubling)', 'tank1.apply_flux(flux_L*A, sign=-1)\ntank2.apply_flux(flux_R*A, sign=+1)'], + ['Phase 6', 'Advance time', 't += dt, step += 1'], + ['Phase 7', 'Record history', 'Store P1, T1, P2, T2, pipe.W'], + ], + col_widths=[2, 4, 9], + ) + + doc.add_page_break() + + # ======================================================================== + # 8. PIPE STEP DETAIL + # ======================================================================== + doc.add_heading('8 Pipe Single-Step Update Detail (pipe.step)', level=1) + + doc.add_heading('8.1 Flux Computation and State Update', level=2) + doc.add_paragraph( + 'The pipe.step(flux_L, flux_R, dt) method is the core of the 1D pipe model. ' + 'It receives the two boundary fluxes from the solver, internally computes N-1 ' + 'interior interface fluxes, and then updates all cells using the finite-volume formula. ' + 'The detailed steps:' + ) + + doc.add_paragraph( + '(1) Snapshot current state: W_snap = W.copy() (ensures same time level)\n' + '\n' + '(2) Compute N-1 interior interface fluxes:\n' + ' For k = 0, 1, ..., N-2:\n' + ' flux_int[:, k] = hll_flux(W_snap[:, k], W_snap[:, k+1], gamma)\n' + '\n' + '(3) Update first cell (Cell 0):\n' + ' Left face = flux_L (from Tank 1 ghost cell)\n' + ' Right face = flux_int[:, 0]\n' + ' W[:, 0] = W_snap[:, 0] - (dt/dx) * (flux_int[:, 0] - flux_L)\n' + '\n' + '(4) Update interior cells (Cell 1 to Cell N-2):\n' + ' W[:, i] = W_snap[:, i] - (dt/dx) * (flux_int[:, i] - flux_int[:, i-1])\n' + '\n' + '(5) Update last cell (Cell N-1):\n' + ' Left face = flux_int[:, N-2]\n' + ' Right face = flux_R (from Tank 2 ghost cell)\n' + ' W[:, N-1] = W_snap[:, N-1] - (dt/dx) * (flux_R - flux_int[:, N-2])' + ) + + doc.add_heading('8.2 Spatial Relationship of Fluxes and Cells', level=2) + add_code_block(doc, + ' flux_L flux_int[0] flux_int[1] flux_int[N-2] flux_R\n' + ' | | | ... | |\n' + ' v v v v v\n' + ' | Cell 0 | Cell 1 | Cell 2 | ... | Cell N-1 |\n' + ' | W[:,0] | W[:,1] | W[:,2] | | W[:,N-1] |\n' + ) + + doc.add_page_break() + + # ======================================================================== + # 9. FRICTION + # ======================================================================== + doc.add_heading('9 Wall Friction: Source Term Method', level=1) + + doc.add_heading('9.1 Modified Governing Equations', level=2) + doc.add_paragraph( + 'With wall friction, the 1D Euler equations become a system with source terms:' + ) + add_equation(doc, 'dW/dt + dF(W)/dx = S') + doc.add_paragraph('Where the source term vector S is:') + add_equation(doc, 'S = [0, -f/D * rho*u*|u|/2, 0]^T') + doc.add_paragraph( + 'Meaning of each component:\n' + '- S[0] = 0: Friction does not create or destroy mass\n' + '- S[1] = -f/D * rho*u*|u|/2: Darcy-Weisbach wall friction force (per unit volume)\n' + ' -- f is the Darcy friction factor (dimensionless)\n' + ' -- D is the pipe inner diameter\n' + ' -- |u| ensures the drag direction always opposes the flow direction\n' + '- S[2] = 0: Under the adiabatic wall assumption, kinetic energy dissipated by friction ' + 'is entirely converted to internal energy; total energy (internal + kinetic) remains unchanged\n' + '\n' + 'Note: S[2] = 0 means wall friction does not change the system total energy. ' + 'Friction decelerates the flow (momentum decreases), but the lost kinetic energy ' + 'is converted to internal energy via frictional heating; their sum remains constant. ' + 'Therefore, even with friction enabled, the system total energy remains strictly conserved.' + ) + + doc.add_heading('9.2 Operator Splitting', level=2) + doc.add_paragraph( + 'To maintain code clarity and modularity, friction source terms are handled via ' + 'operator splitting: within each time step, first complete the source-free flux update ' + '(Euler equation part), then separately apply the source term. ' + 'This is equivalent to Lie Splitting:' + ) + add_equation(doc, 'W* = W^(n) - (dt/dx)*[F_{i+1/2} - F_{i-1/2}] (flux step)') + add_equation(doc, 'W^(n+1)[1] = W*[1] + dt * S[1] (source step, momentum only)') + doc.add_paragraph( + 'Since S[0] = S[2] = 0, the source step only updates the momentum component W[1] = rho*u. ' + 'Implementation:' + ) + add_code_block(doc, + 'if self.mu > 0:\n' + ' rho_s = self.W[0, :]\n' + ' u_s = self.W[1, :] / rho_s\n' + ' abs_u = np.abs(u_s)\n' + ' Re = rho_s * abs_u * self.D / self.mu\n' + ' f = darcy_friction_factor(Re, self.eps_D)\n' + ' S_mom = -f / self.D * rho_s * u_s * abs_u / 2.0\n' + ' self.W[1, :] += dt * S_mom') + + doc.add_heading('9.3 Darcy Friction Factor Calculation', level=2) + doc.add_paragraph( + 'The friction factor f depends on the Reynolds number Re and relative wall roughness eps/D:' + ) + add_equation(doc, 'Re = rho * |u| * D / mu') + doc.add_paragraph( + 'Where mu is the dynamic viscosity. Depending on Re, three regimes are used:' + ) + + make_table(doc, + ['Flow Regime', 'Re Range', 'Friction Factor Formula'], + [ + ['Laminar', 'Re < 2300', 'f = 64 / Re'], + ['Transition', '2300 <= Re <= 4000', 'Linear blend of laminar and turbulent:\n' + 'f = (1-alpha)*f_lam + alpha*f_turb\nalpha = (Re - 2300) / 1700'], + ['Turbulent', 'Re > 4000', 'Colebrook-White implicit equation:\n' + '1/sqrt(f) = -2*log10(eps/(3.7*D) + 2.51/(Re*sqrt(f)))\n' + 'Solved via Swamee-Jain initial guess + 10 fixed-point iterations'], + ], + col_widths=[2.5, 4, 9], + ) + + doc.add_paragraph( + 'The current default configuration uses smooth pipe walls (ROUGHNESS = 0), ' + 'in which case the Colebrook-White equation simplifies to the smooth-pipe implicit friction law.' + ) + + doc.add_page_break() + + # ======================================================================== + # 10. CONSERVATION PROOF + # ======================================================================== + doc.add_heading('10 Conservation Analysis and Verification', level=1) + + doc.add_heading('10.1 Total System Mass', level=2) + add_equation(doc, 'M_total = m_tank1 + m_tank2 + SUM_i(rho_i * A * dx)') + doc.add_paragraph( + 'Where the summation runs over all N pipe cells. Due to the flux doubling mechanism, ' + 'within each time step:\n' + '- Tank 1 mass loss = flux_L[0] * A * dt\n' + '- Cell 0 mass gain = (flux_L[0] - flux_int[0]) * (A*dt/dx) * dx = difference\n' + '- ... internal cell flux differences exactly "pass" mass through ...\n' + '- Cell N-1 mass gain = (flux_int[N-2] - flux_R[0]) * ...\n' + '- Tank 2 mass gain = flux_R[0] * A * dt\n' + '\n' + 'This is a Telescoping Sum: all intermediate terms cancel, and the boundary terms ' + 'are exactly offset by the tank updates. Therefore M_total is strictly invariant.' + ) + + doc.add_heading('10.2 Total System Energy', level=2) + add_equation(doc, 'E_total = U_tank1 + U_tank2 + SUM_i((rho*E)_i * A * dx)') + doc.add_paragraph( + 'The analysis is entirely analogous to mass conservation. ' + 'Note: friction source term S[2] = 0, so friction does not affect energy conservation.' + ) + + doc.add_heading('10.3 Numerical Verification Results', level=2) + doc.add_paragraph( + 'Conservation checks are performed in both main.py and tests/test_integration.py:' + ) + make_table(doc, + ['Conserved Quantity', 'Initial Value', 'Final Value', 'Relative Error'], + [ + ['Total mass', '5.226485 x 10^2 kg', '5.226485 x 10^2 kg', '~6.5 x 10^(-16)'], + ['Total energy', '1.125001 x 10^8 J', '1.125001 x 10^8 J', '~1.1 x 10^(-15)'], + ], + ) + doc.add_paragraph( + 'Relative errors are at the 10^(-15) level, i.e., IEEE 754 double-precision ' + 'machine epsilon, confirming the correctness of the flux doubling mechanism.' + ) + + doc.add_page_break() + + # ======================================================================== + # 11. SUMMARY OF DATA FLOW + # ======================================================================== + doc.add_heading('11 Data Flow Summary', level=1) + + doc.add_paragraph( + 'The following summarizes the data flow between modules in one complete time step:' + ) + add_code_block(doc, + '+-----------------------------------------------------------+\n' + '| solver.run() |\n' + '| |\n' + '| (1) CFL: a_max = pipe.max_wave_speed() |\n' + '| dt = CFL * dx / a_max |\n' + '| |\n' + '| (2) Ghost: W_gL = tank1.ghost_state() |\n' + '| W_gR = tank2.ghost_state() |\n' + '| |\n' + '| (3) Boundary flux: |\n' + '| fL = hll_flux(W_gL, pipe.W[:,0]) <-- riemann.py |\n' + '| fR = hll_flux(pipe.W[:,-1], W_gR) <-- riemann.py |\n' + '| |\n' + '| (4) Pipe update: |\n' + '| pipe.step(fL, fR, dt) |\n' + '| +-- internal HLL fluxes (riemann.py) |\n' + '| +-- finite-volume update W |\n' + '| +-- friction source term (friction.py) |\n' + '| |\n' + '| (5) Tank update (same fL, fR): |\n' + '| tank1.apply_flux(fL*A, sign=-1) |\n' + '| tank2.apply_flux(fR*A, sign=+1) |\n' + '| |\n' + '| (6)(7) t += dt, record history |\n' + '+-----------------------------------------------------------+\n' + ) + + doc.add_page_break() + + # ======================================================================== + # 12. PARAMETERS + # ======================================================================== + doc.add_heading('12 Simulation Parameter Summary', level=1) + + doc.add_heading('12.1 Gas Properties', level=2) + make_table(doc, + ['Parameter', 'Symbol', 'Value', 'Unit'], + [ + ['Heat capacity ratio', 'gamma', '1.4', '--'], + ['Gas constant', 'R', '287.0', 'J/(kg*K)'], + ['Dynamic viscosity', 'mu', '1.8 x 10^(-5)', 'Pa*s'], + ], + ) + + doc.add_heading('12.2 Geometry and Initial Conditions', level=2) + make_table(doc, + ['Parameter', 'Symbol', 'Value', 'Unit'], + [ + ['High-pressure tank volume', 'V1', '5.0', 'm^3'], + ['High-pressure tank initial pressure', 'P1', '5.0', 'MPa'], + ['High-pressure tank initial temperature', 'T1', '300.0', 'K'], + ['Low-pressure tank volume', 'V2', '10.0', 'm^3'], + ['Low-pressure tank initial pressure', 'P2', '2.0', 'MPa'], + ['Low-pressure tank initial temperature', 'T2', '300.0', 'K'], + ['Pipe length', 'L', '1.0', 'm'], + ['Pipe inner diameter', 'D', '5.0', 'mm'], + ['Wall roughness', 'eps', '0.0', 'm (smooth pipe)'], + ], + ) + + doc.add_heading('12.3 Numerical Control', level=2) + make_table(doc, + ['Parameter', 'Symbol', 'Value', 'Description'], + [ + ['Number of cells', 'N', '20', 'Pipe uniformly divided into 20 finite-volume cells'], + ['Cell size', 'dx', '50 mm', '= L / N'], + ['CFL number', 'CFL', '0.5', 'Time-step safety factor'], + ['Simulation end time', 't_end', '0.1', 'seconds'], + ], + ) + + # ======================================================================== + # SAVE + # ======================================================================== + out_path = os.path.join('docs', 'pipe_system_simulation_technical_doc.docx') + os.makedirs('docs', exist_ok=True) + doc.save(out_path) + print(f'Document saved to: {out_path}') + return out_path + + +if __name__ == '__main__': + build_document() diff --git a/results/.gitkeep b/results/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 0000000..aa0f361 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,3 @@ +# scripts + +Automation, preprocessing, and postprocessing scripts. diff --git a/src/README.md b/src/README.md new file mode 100644 index 0000000..f09c834 --- /dev/null +++ b/src/README.md @@ -0,0 +1,3 @@ +# src + +Source code for pipe system simulation models and utilities. diff --git a/src/config.py b/src/config.py new file mode 100644 index 0000000..7ed3dae --- /dev/null +++ b/src/config.py @@ -0,0 +1,52 @@ +""" +Physical, geometric, and numerical constants for the 0D-1D tank-pipe +blowdown MVP. Pure data module — no functions, no side effects. +""" + +# ---------- Gas properties (ideal air-like) ---------- +GAMMA = 1.4 +R_GAS = 287.0 # J / (kg K) + +# ---------- High-pressure tank (upstream, Tank 1) ---------- +V1 = 5.0 # m^3 +P1_INIT = 10e6 # Pa (10 MPa) +T1_INIT = 300.0 # K + +# ---------- Low-pressure tank (downstream, Tank 2) ---------- +V2 = 10.0 # m^3 +P2_INIT = 2e6 # Pa (2 MPa) +T2_INIT = 300.0 # K + +# ---------- Pipe geometry ---------- +L = 1.0 # m +D = 5e-3 # m (5 mm) +N_CELLS = 20 # number of finite-volume cells + +# ---------- Friction ---------- +MU = 1.8e-5 # Pa·s (dynamic viscosity of air at ~300 K) +ROUGHNESS = 0.0 # m (absolute wall roughness; 0 = smooth pipe) + +# ---------- Simulation control ---------- +T_END = 0.1 # s +CFL = 0.5 +RIEMANN_SOLVER = "roe" # "hll" or "roe" + +# ---------- Output & animation ---------- +ANIMATION_STRIDE = 10 # keep every Nth frame in the GIF +OUTPUT_DIR = "results" + +# ---------- Parameter validation (per spec §6.1) ---------- +assert GAMMA > 1, "GAMMA must be > 1" +assert R_GAS > 0, "R_GAS must be > 0" +assert V1 > 0 and V2 > 0, "tank volumes must be > 0" +assert L > 0, "L must be > 0" +assert D > 0, "D must be > 0" +assert N_CELLS >= 2, "N_CELLS must be >= 2" +assert P1_INIT > 0 and P2_INIT > 0, "initial pressures must be > 0" +assert T1_INIT > 0 and T2_INIT > 0, "initial temperatures must be > 0" +assert MU >= 0, "MU must be >= 0" +assert ROUGHNESS >= 0, "ROUGHNESS must be >= 0" +assert 0 < CFL <= 1, "CFL must be in (0, 1]" +assert RIEMANN_SOLVER in ("hll", "roe"), "RIEMANN_SOLVER must be 'hll' or 'roe'" +assert T_END > 0, "T_END must be > 0" +assert ANIMATION_STRIDE >= 1, "ANIMATION_STRIDE must be >= 1" diff --git a/src/cryo_tank/__init__.py b/src/cryo_tank/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/cryo_tank/config.py b/src/cryo_tank/config.py new file mode 100644 index 0000000..01f0b71 --- /dev/null +++ b/src/cryo_tank/config.py @@ -0,0 +1,63 @@ +# 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 diff --git a/src/cryo_tank/heat_leak.py b/src/cryo_tank/heat_leak.py new file mode 100644 index 0000000..6d1ff62 --- /dev/null +++ b/src/cryo_tank/heat_leak.py @@ -0,0 +1,70 @@ +# 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 diff --git a/src/cryo_tank/main.py b/src/cryo_tank/main.py new file mode 100644 index 0000000..c64e31e --- /dev/null +++ b/src/cryo_tank/main.py @@ -0,0 +1,80 @@ +# 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, plot_pressure, +) + + +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")) + plot_pressure(history, os.path.join(OUTPUT_DIR, "cryo_tank_pressure.png")) + + print(f"\nOutputs written to {OUTPUT_DIR}/") + + +if __name__ == "__main__": + main() diff --git a/src/cryo_tank/output.py b/src/cryo_tank/output.py new file mode 100644 index 0000000..58030a1 --- /dev/null +++ b/src/cryo_tank/output.py @@ -0,0 +1,112 @@ +# 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) + + +def plot_pressure(history, path): + """Plot tank pressure (P_total, P_N2, P_He) vs time.""" + fig, ax = plt.subplots(figsize=(10, 5)) + t = history['t'] + + ax.plot(t, history['P_total'] / 1e6, label='P_total', linewidth=2) + ax.plot(t, history['P_N2'] / 1e6, label='P_N2', linestyle='--') + ax.plot(t, history['P_He'] / 1e6, label='P_He', linestyle='--') + ax.set_xlabel('Time [s]') + ax.set_ylabel('Pressure [MPa]') + ax.set_title('Tank pressure vs time') + ax.grid(True) + ax.legend() + fig.tight_layout() + fig.savefig(path, dpi=120) + plt.close(fig) diff --git a/src/cryo_tank/properties.py b/src/cryo_tank/properties.py new file mode 100644 index 0000000..8d23d0d --- /dev/null +++ b/src/cryo_tank/properties.py @@ -0,0 +1,121 @@ +# 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 diff --git a/src/cryo_tank/solver.py b/src/cryo_tank/solver.py new file mode 100644 index 0000000..90fa3a3 --- /dev/null +++ b/src/cryo_tank/solver.py @@ -0,0 +1,130 @@ +# 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), + 'P_N2': np.zeros(n), + 'P_He': np.zeros(n), + 'P_total': 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'] + history['P_N2'][i] = info['P_N2'] + history['P_He'][i] = info['P_He'] + history['P_total'][i] = info['P_N2'] + info['P_He'] + + # 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 diff --git a/src/cryo_tank/tank_model.py b/src/cryo_tank/tank_model.py new file mode 100644 index 0000000..95c311b --- /dev/null +++ b/src/cryo_tank/tank_model.py @@ -0,0 +1,282 @@ +# 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 + # Use ideal gas law for N2 mass (consistent with derive/rhs which + # treat ullage N2 as ideal gas: P_N2 = m_N2 * R_N2 * T / V) + P_N2_0 = prop.n2_sat_pressure(T_init) + self.m_N2_ull = P_N2_0 * V_ull_0 / (R_N2 * T_init) # 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 + + # Solve T_ull from ullage 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 + 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] + + 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 diff --git a/src/friction.py b/src/friction.py new file mode 100644 index 0000000..20c227f --- /dev/null +++ b/src/friction.py @@ -0,0 +1,74 @@ +# src/friction.py +""" +Darcy-Weisbach friction factor calculation. + +Supports: + - Laminar: f = 64 / Re (Re < 2300) + - Turbulent: Colebrook-White implicit equation (Re > 4000) + - Transition: linear blend between laminar & turbulent (2300 <= Re <= 4000) +""" + +import numpy as np + + +def _colebrook_white(Re, eps_D, n_iter=10): + """ + Solve the Colebrook-White equation for Darcy friction factor f: + 1/sqrt(f) = -2 log10( eps_D/3.7 + 2.51/(Re*sqrt(f)) ) + + Uses fixed-point iteration seeded with the Swamee-Jain approximation. + """ + # Swamee-Jain initial guess (explicit approximation) + A = eps_D / 3.7 + B = 2.51 / Re + f = 0.25 / (np.log10(A + B / np.sqrt(0.02))) ** 2 + + for _ in range(n_iter): + f = 0.25 / (np.log10(A + B / np.sqrt(f))) ** 2 + + return f + + +def darcy_friction_factor(Re, eps_D): + """ + Compute Darcy-Weisbach friction factor for a given Reynolds number + and relative roughness eps/D. + + Parameters + ---------- + Re : float or ndarray + Reynolds number (ρ|u|D/μ). Values <= 0 return 0 (no flow). + eps_D : float + Relative roughness ε/D (dimensionless). + + Returns + ------- + f : same shape as Re + Darcy friction factor. + """ + Re = np.asarray(Re, dtype=float) + scalar = Re.ndim == 0 + Re = np.atleast_1d(Re) + + f = np.zeros_like(Re) + + lam = Re < 2300 + turb = Re > 4000 + trans = ~lam & ~turb # 2300 <= Re <= 4000 + + # Laminar: f = 64/Re (avoid division by zero for Re~0) + Re_lam = np.where(Re > 1e-12, Re, 1e-12) + f[lam] = 64.0 / Re_lam[lam] + + # Turbulent: Colebrook-White + if np.any(turb): + f[turb] = _colebrook_white(Re[turb], eps_D) + + # Transition: linear blend + if np.any(trans): + f_lam = 64.0 / Re_lam[trans] + f_turb = _colebrook_white(Re[trans], eps_D) + alpha = (Re[trans] - 2300.0) / 1700.0 # 0 at Re=2300, 1 at Re=4000 + f[trans] = (1.0 - alpha) * f_lam + alpha * f_turb + + return float(f[0]) if scalar else f diff --git a/src/main.py b/src/main.py new file mode 100644 index 0000000..5772503 --- /dev/null +++ b/src/main.py @@ -0,0 +1,131 @@ +# src/main.py +""" +Entry point: assemble tanks + pipe from config constants, run the solver, +verify total mass/energy conservation, persist history, and generate +plots + animation. + +Run from project root: + python3 src/main.py +""" +import os +import sys + +# Ensure imports work when running from project root +_HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, _HERE) + +import numpy as np + +from config import ( + GAMMA, R_GAS, + V1, P1_INIT, T1_INIT, + V2, P2_INIT, T2_INIT, + L, D, N_CELLS, + MU, ROUGHNESS, + T_END, CFL, RIEMANN_SOLVER, + ANIMATION_STRIDE, OUTPUT_DIR, +) +from tank import Tank +from pipe import Pipe +from solver import run +from output import ( + save_history, + plot_tank_pressure, + plot_tank_temperature, + plot_pipe_final_profiles, + make_pipe_animation, + write_summary_report, +) + + +def _total_mass(tank1, tank2, pipe): + pipe_mass = float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx)) + return tank1.mass + tank2.mass + pipe_mass + + +def _total_energy(tank1, tank2, pipe): + pipe_energy = float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx)) + return tank1.U + tank2.U + pipe_energy + + +def main(): + os.makedirs(OUTPUT_DIR, exist_ok=True) + + # --- Assemble --- + tank1 = Tank(V=V1, P_init=P1_INIT, T_init=T1_INIT, gamma=GAMMA, R_gas=R_GAS) + tank2 = Tank(V=V2, P_init=P2_INIT, T_init=T2_INIT, gamma=GAMMA, R_gas=R_GAS) + pipe = Pipe(L=L, D=D, N=N_CELLS, P_init=P2_INIT, T_init=T2_INIT, + gamma=GAMMA, R_gas=R_GAS, mu=MU, roughness=ROUGHNESS, + riemann_solver=RIEMANN_SOLVER) + + m_init = _total_mass(tank1, tank2, pipe) + U_init = _total_energy(tank1, tank2, pipe) + print(f"Initial total mass: {m_init:.6e} kg") + print(f"Initial total energy: {U_init:.6e} J") + print(f"Initial P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa") + print(f"Pipe: L={L} m, D={D*1e3:.1f} mm, N={N_CELLS} cells, dx={pipe.dx*1e3:.1f} mm") + print(f"Riemann solver: {RIEMANN_SOLVER.upper()}") + if MU > 0: + print(f"Friction: mu={MU:.2e} Pa·s, roughness={ROUGHNESS:.2e} m (eps/D={ROUGHNESS/D:.4f})") + else: + print("Friction: OFF") + print(f"Running to t_end={T_END} s with CFL={CFL}...") + print() + + # --- Run --- + history = run(tank1, tank2, pipe, + t_end=T_END, cfl=CFL, + verbose=True, log_every=200) + + n_steps = len(history['t']) + print() + print(f"Simulation complete: {n_steps} steps") + + # --- Conservation sanity check (per spec §6.1, §8) --- + m_final = _total_mass(tank1, tank2, pipe) + U_final = _total_energy(tank1, tank2, pipe) + rel_err_m = abs(m_final - m_init) / m_init + rel_err_U = abs(U_final - U_init) / U_init + print(f"Final total mass: {m_final:.6e} kg (rel err = {rel_err_m:.2e})") + print(f"Final total energy: {U_final:.6e} J (rel err = {rel_err_U:.2e})") + print(f"Final P1 = {tank1.P/1e6:.3f} MPa, P2 = {tank2.P/1e6:.3f} MPa") + assert rel_err_m < 1e-10, f"Total mass not conserved: rel_err={rel_err_m:.2e}" + assert rel_err_U < 1e-10, f"Total energy not conserved: rel_err={rel_err_U:.2e}" + + # --- Persist + visualize --- + save_history(history, pipe, + os.path.join(OUTPUT_DIR, "history.npz"), + GAMMA, R_GAS) + plot_tank_pressure(history, + os.path.join(OUTPUT_DIR, "tank_pressure.png")) + plot_tank_temperature(history, + os.path.join(OUTPUT_DIR, "tank_temperature.png")) + plot_pipe_final_profiles(history, pipe, + os.path.join(OUTPUT_DIR, "pipe_final_profiles.png"), + GAMMA, R_GAS) + make_pipe_animation(history, pipe, + os.path.join(OUTPUT_DIR, "pipe_animation.gif"), + GAMMA, R_GAS, stride=ANIMATION_STRIDE) + write_summary_report( + history, pipe, + os.path.join(OUTPUT_DIR, "summary_report.html"), + GAMMA, R_GAS, + config={ + 'V1': V1, 'P1_INIT': P1_INIT, 'T1_INIT': T1_INIT, + 'V2': V2, 'P2_INIT': P2_INIT, 'T2_INIT': T2_INIT, + 'L': L, 'D': D, 'N_CELLS': N_CELLS, + 'T_END': T_END, 'CFL': CFL, + }, + ) + + print(f"Outputs written to {OUTPUT_DIR}/") + print(f" - history.npz") + print(f" - tank_pressure.png") + print(f" - tank_temperature.png") + print(f" - pipe_final_profiles.png") + print(f" - pipe_animation.gif") + print(f" - summary_report.html") + + +if __name__ == "__main__": + main() diff --git a/src/output.py b/src/output.py new file mode 100644 index 0000000..99a8ab0 --- /dev/null +++ b/src/output.py @@ -0,0 +1,300 @@ +# src/output.py +""" +Output helpers: persistence (.npz), static plots (.png/.html), animation (.gif). +Uses matplotlib's Agg backend so it works in headless environments. +The PillowWriter is used for GIF output to avoid an ffmpeg dependency. +""" +import html +import os + +import numpy as np +import matplotlib +matplotlib.use("Agg") # headless-safe; must be set before pyplot import +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation, PillowWriter + + +def _pipe_primitives(history, gamma, R_gas): + W_hist = history['W_hist'] + rho = W_hist[:, 0, :] + u = W_hist[:, 1, :] / rho + P = (gamma - 1) * (W_hist[:, 2, :] - 0.5 * rho * u ** 2) + T = P / (rho * R_gas) + a = np.sqrt(gamma * P / rho) + Ma = u / a + return rho, u, P, T, a, Ma + + +def save_history(history, pipe, path, gamma, R_gas): + """ + Persist the full simulation history + pipe geometry to a .npz file. + + Loadable later with: + d = np.load("results/history.npz") + rho = d['W_hist'][:, 0, :] + u = d['W_hist'][:, 1, :] / rho + P = (d['gamma'] - 1) * (d['W_hist'][:, 2, :] - 0.5 * rho * u**2) + """ + dirname = os.path.dirname(path) + if dirname: + os.makedirs(dirname, exist_ok=True) + np.savez_compressed( + path, + t=history['t'], + P1=history['P1'], T1=history['T1'], + P2=history['P2'], T2=history['T2'], + W_hist=history['W_hist'], + x=pipe.x_centers, + dx=pipe.dx, + area=pipe.area, + gamma=gamma, + R_gas=R_gas, + ) + + +def plot_tank_pressure(history, path): + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(history['t'], history['P1'] / 1e6, label="Tank 1 (high pressure)") + ax.plot(history['t'], history['P2'] / 1e6, label="Tank 2 (low pressure)") + ax.set_xlabel("Time [s]") + ax.set_ylabel("Pressure [MPa]") + ax.set_title("Tank pressures vs time") + ax.grid(True) + ax.legend() + fig.tight_layout() + fig.savefig(path, dpi=120) + plt.close(fig) + + +def plot_tank_temperature(history, path): + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(history['t'], history['T1'], label="Tank 1 (high pressure)") + ax.plot(history['t'], history['T2'], label="Tank 2 (low pressure)") + 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 make_pipe_animation(history, pipe, path, gamma, R_gas, stride=10): + """ + Render a GIF of the pipe's P(x), u(x), T(x) evolution over time. + Uses PillowWriter so no ffmpeg is needed. + """ + t = history['t'] + x = pipe.x_centers + n_steps = history['W_hist'].shape[0] + + # Frame indices: every `stride`th snapshot, plus the final one + frames = list(range(0, n_steps, stride)) + if frames[-1] != n_steps - 1: + frames.append(n_steps - 1) + + # Precompute primitives for all frames in one vectorized pass + _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas) + + fig, axes = plt.subplots(4, 1, figsize=(10, 12), sharex=True) + + # Pressure subplot + line_P, = axes[0].plot(x, P[0] / 1e6) + axes[0].set_ylabel("P [MPa]") + axes[0].set_ylim(P.min() / 1e6 * 0.95, P.max() / 1e6 * 1.05) + axes[0].grid(True) + + # Velocity subplot + line_u, = axes[1].plot(x, u[0]) + axes[1].set_ylabel("u [m/s]") + u_min, u_max = float(u.min()), float(u.max()) + pad = max(1.0, 0.05 * (u_max - u_min) if u_max > u_min else 1.0) + axes[1].set_ylim(u_min - pad, u_max + pad) + axes[1].grid(True) + + # Temperature subplot + line_T, = axes[2].plot(x, T[0]) + axes[2].set_ylabel("T [K]") + axes[2].set_ylim(T.min() * 0.95, T.max() * 1.05) + axes[2].grid(True) + + # Mach number subplot + line_Ma, = axes[3].plot(x, Ma[0]) + axes[3].set_ylabel("Mach [-]") + axes[3].set_xlabel("x [m]") + Ma_min, Ma_max = float(Ma.min()), float(Ma.max()) + pad_Ma = max(0.05, 0.05 * (Ma_max - Ma_min) if Ma_max > Ma_min else 0.05) + axes[3].set_ylim(Ma_min - pad_Ma, Ma_max + pad_Ma) + axes[3].grid(True) + + title = fig.suptitle("") + + def update(frame_idx): + line_P.set_ydata(P[frame_idx] / 1e6) + line_u.set_ydata(u[frame_idx]) + line_T.set_ydata(T[frame_idx]) + line_Ma.set_ydata(Ma[frame_idx]) + title.set_text( + f"t = {t[frame_idx]:.5f} s (step {frame_idx + 1}/{n_steps})" + ) + return line_P, line_u, line_T, line_Ma, title + + anim = FuncAnimation(fig, update, frames=frames, interval=50, blit=False) + + dirname = os.path.dirname(path) + if dirname: + os.makedirs(dirname, exist_ok=True) + anim.save(path, writer=PillowWriter(fps=20)) + plt.close(fig) + + +def plot_pipe_final_profiles(history, pipe, path, gamma, R_gas): + """Plot final pipe P(x), u(x), T(x), and Mach(x) as a static PNG.""" + x = pipe.x_centers + _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas) + + fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True) + axes = axes.ravel() + + axes[0].plot(x, P[-1] / 1e6) + axes[0].set_ylabel("P [MPa]") + axes[0].set_title("Final pressure profile") + axes[0].grid(True) + + axes[1].plot(x, u[-1]) + axes[1].set_ylabel("u [m/s]") + axes[1].set_title("Final velocity profile") + axes[1].grid(True) + + axes[2].plot(x, T[-1]) + axes[2].set_xlabel("x [m]") + axes[2].set_ylabel("T [K]") + axes[2].set_title("Final temperature profile") + axes[2].grid(True) + + axes[3].plot(x, Ma[-1]) + axes[3].axhline(1.0, color="r", linestyle="--", linewidth=1, label="Mach 1") + axes[3].set_xlabel("x [m]") + axes[3].set_ylabel("Mach [-]") + axes[3].set_title("Final Mach profile") + axes[3].grid(True) + axes[3].legend() + + fig.tight_layout() + fig.savefig(path, dpi=120) + plt.close(fig) + + +def write_summary_report(history, pipe, path, gamma, R_gas, config): + """Write an HTML summary report for the latest simulation outputs.""" + t = history['t'] + P1 = history['P1'] + T1 = history['T1'] + P2 = history['P2'] + T2 = history['T2'] + W_hist = history['W_hist'] + x = pipe.x_centers + + _, u, P, T, _, Ma = _pipe_primitives(history, gamma, R_gas) + + V1 = config['V1'] + V2 = config['V2'] + + m1_init = P1[0] * V1 / (R_gas * T1[0]) + m1_final = P1[-1] * V1 / (R_gas * T1[-1]) + m2_init = P2[0] * V2 / (R_gas * T2[0]) + m2_final = P2[-1] * V2 / (R_gas * T2[-1]) + mpipe_init = float(np.sum(W_hist[0, 0, :] * pipe.area * pipe.dx)) + mpipe_final = float(np.sum(W_hist[-1, 0, :] * pipe.area * pipe.dx)) + m_total_init = m1_init + m2_init + mpipe_init + m_total_final = m1_final + m2_final + mpipe_final + rel_err_m = abs(m_total_final - m_total_init) / m_total_init + + U1_init = P1[0] * V1 / (gamma - 1) + U1_final = P1[-1] * V1 / (gamma - 1) + U2_init = P2[0] * V2 / (gamma - 1) + U2_final = P2[-1] * V2 / (gamma - 1) + Upipe_init = float(np.sum(W_hist[0, 2, :] * pipe.area * pipe.dx)) + Upipe_final = float(np.sum(W_hist[-1, 2, :] * pipe.area * pipe.dx)) + U_total_init = U1_init + U2_init + Upipe_init + U_total_final = U1_final + U2_final + Upipe_final + rel_err_U = abs(U_total_final - U_total_init) / U_total_init + + idx_ma = np.unravel_index(np.argmax(Ma), Ma.shape) + idx_p = np.unravel_index(np.argmax(P), P.shape) + + sections = [ + ("Simulation setup", [ + ("Gamma", f"{gamma:.3f}"), + ("R_gas", f"{R_gas:.3f} J/(kg K)"), + ("Tank 1", f"V={V1:.3f} m^3, P0={config['P1_INIT']/1e6:.6f} MPa, T0={config['T1_INIT']:.3f} K"), + ("Tank 2", f"V={V2:.3f} m^3, P0={config['P2_INIT']/1e6:.6f} MPa, T0={config['T2_INIT']:.3f} K"), + ("Pipe", f"L={config['L']:.3f} m, D={config['D']*1e3:.3f} mm, N={config['N_CELLS']}"), + ("Run control", f"t_end={config['T_END']:.6f} s, CFL={config['CFL']:.3f}, steps={len(t)}"), + ]), + ("Tank states", [ + ("Tank 1 pressure", f"{P1[0]/1e6:.6f} -> {P1[-1]/1e6:.6f} MPa"), + ("Tank 2 pressure", f"{P2[0]/1e6:.6f} -> {P2[-1]/1e6:.6f} MPa"), + ("Tank 1 temperature", f"{T1[0]:.6f} -> {T1[-1]:.6f} K"), + ("Tank 2 temperature", f"{T2[0]:.6f} -> {T2[-1]:.6f} K"), + ]), + ("Pipe extrema", [ + ("Max pressure", f"{P[idx_p]/1e6:.6f} MPa at t={t[idx_p[0]]:.6e} s, x={x[idx_p[1]]:.6f} m"), + ("Max velocity", f"{u.max():.6f} m/s"), + ("Min / max temperature", f"{T.min():.6f} / {T.max():.6f} K"), + ("Max Mach", f"{Ma[idx_ma]:.6f} at t={t[idx_ma[0]]:.6e} s, x={x[idx_ma[1]]:.6f} m"), + ("Final Mach range", f"{Ma[-1].min():.6f} -> {Ma[-1].max():.6f}"), + ("Final supersonic cells", f"{int(np.sum(Ma[-1] > 1.0))} / {Ma.shape[1]}"), + ]), + ("Conservation check", [ + ("Total mass", f"{m_total_init:.12e} -> {m_total_final:.12e} kg (rel err {rel_err_m:.3e})"), + ("Total energy", f"{U_total_init:.12e} -> {U_total_final:.12e} J (rel err {rel_err_U:.3e})"), + ]), + ] + + parts = [ + "", + "", + "", + "", + "Pipe system simulation summary", + "", + "", + "", + "

Pipe system simulation summary

", + f"

Generated from results/history.npz. Final simulation time: {t[-1]:.6f} s.

", + ] + + for title, rows in sections: + parts.append(f"

{html.escape(title)}

") + parts.append("") + for key, value in rows: + parts.append( + f"" + ) + parts.append("
{html.escape(str(key))}{html.escape(str(value))}
") + + parts.extend([ + "

Figures

", + "

Tank pressure history

", + "

Tank temperature history

", + "

Final pipe profiles

", + "

Animation: pipe_animation.gif

", + "", + "", + ]) + + dirname = os.path.dirname(path) + if dirname: + os.makedirs(dirname, exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + f.write("\n".join(parts)) diff --git a/src/pipe.py b/src/pipe.py new file mode 100644 index 0000000..1ef0c36 --- /dev/null +++ b/src/pipe.py @@ -0,0 +1,126 @@ +# src/pipe.py +""" +1D finite-volume pipe for compressible Euler equations: + dW/dt + dF(W)/dx = 0 + W = [rho, rho*u, rho*E], F = [rho*u, rho*u^2 + P, u*(rho*E + P)] + +Discretization: + - N uniform cells, cell-averaged piecewise-constant reconstruction + - HLL numerical flux at all interior interfaces + - Boundary (tank-side) interface fluxes are provided by the caller + via step(flux_L, flux_R, dt) +""" + +import numpy as np +from riemann import hll_flux, get_riemann_solver +from friction import darcy_friction_factor + + +class Pipe: + def __init__(self, L, D, N, P_init, T_init, gamma, R_gas, + mu=0.0, roughness=0.0, riemann_solver="hll"): + self.L = L + self.D = D + self.N = N + self.dx = L / N + self.area = np.pi * (D / 2) ** 2 + self.gamma = gamma + self.R = R_gas + self.mu = mu # dynamic viscosity [Pa·s] + self.roughness = roughness # absolute wall roughness [m] + self.eps_D = roughness / D if D > 0 else 0.0 # relative roughness + self._flux_fn = get_riemann_solver(riemann_solver) + self.x_centers = np.linspace(self.dx / 2, L - self.dx / 2, N) + + # Uniform initial state, u = 0 + rho = P_init / (R_gas * T_init) + E_density = P_init / (gamma - 1) # since u=0, total energy density = internal + + self.W = np.zeros((3, N)) + self.W[0, :] = rho + self.W[1, :] = 0.0 + self.W[2, :] = E_density + + def primitives(self): + """ + Return (rho, u, P, a) each of shape (N,), computed from W. + """ + rho = self.W[0, :] + u = self.W[1, :] / rho + P = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho * u ** 2) + a = np.sqrt(self.gamma * P / rho) + return rho, u, P, a + + def max_wave_speed(self): + """ + Return max over cells of |u| + a, used for CFL dt calculation. + """ + _, u, _, a = self.primitives() + return float(np.max(np.abs(u) + a)) + + def step(self, flux_L, flux_R, dt): + """ + Advance W by one explicit Euler step. The caller provides the + two boundary interface fluxes (with ghost states already folded + in); internal interface fluxes are computed here with HLL. + + Parameters + ---------- + flux_L, flux_R : np.ndarray of shape (3,) + Numerical fluxes at the leftmost and rightmost interfaces + (cell -1/2 and cell N-1/2, i.e. the tank-facing boundaries). + dt : float + Time-step size. + + Raises + ------ + RuntimeError + If the updated state has any non-positive density or pressure. + """ + N = self.N + W_snap = self.W.copy() + + # Internal fluxes: flux_int[:, k] is the flux at the interface + # between cell k and cell k+1, for k = 0 .. N-2 (total N-1 of them) + flux_int = np.zeros((3, N - 1)) + for k in range(N - 1): + flux_int[:, k] = self._flux_fn(W_snap[:, k], W_snap[:, k + 1], self.gamma) + + # First cell: left face = flux_L, right face = flux_int[:, 0] + self.W[:, 0] = W_snap[:, 0] - (dt / self.dx) * (flux_int[:, 0] - flux_L) + + # Interior cells: left face = flux_int[:, i-1], right face = flux_int[:, i] + for i in range(1, N - 1): + self.W[:, i] = W_snap[:, i] - (dt / self.dx) * (flux_int[:, i] - flux_int[:, i - 1]) + + # Last cell: left face = flux_int[:, N-2], right face = flux_R + self.W[:, N - 1] = W_snap[:, N - 1] - (dt / self.dx) * (flux_R - flux_int[:, N - 2]) + + # --- Friction source term (operator splitting, explicit Euler) --- + # S = [0, -f/D * rho*u*|u|/2, 0] + # Energy source = 0 for adiabatic wall (KE dissipated → internal energy) + if self.mu > 0: + rho_s = self.W[0, :] + u_s = self.W[1, :] / rho_s + abs_u = np.abs(u_s) + Re = rho_s * abs_u * self.D / self.mu + f = darcy_friction_factor(Re, self.eps_D) + S_mom = -f / self.D * rho_s * u_s * abs_u / 2.0 + self.W[1, :] += dt * S_mom + + # Physical-state sanity check + rho_new = self.W[0, :] + if np.any(rho_new <= 0): + bad = np.where(rho_new <= 0)[0] + raise RuntimeError( + f"Non-positive density after pipe step at cells {bad.tolist()}: " + f"rho={rho_new[bad].tolist()}" + ) + u_new = self.W[1, :] / rho_new + P_new = (self.gamma - 1) * (self.W[2, :] - 0.5 * rho_new * u_new ** 2) + if np.any(P_new <= 0): + bad = np.where(P_new <= 0)[0] + raise RuntimeError( + f"Non-positive pressure after pipe step at cells {bad.tolist()}: " + f"P={P_new[bad].tolist()}" + ) diff --git a/src/riemann.py b/src/riemann.py new file mode 100644 index 0000000..0e8a3dc --- /dev/null +++ b/src/riemann.py @@ -0,0 +1,206 @@ +# src/riemann.py +""" +Riemann flux solvers for the 1D compressible Euler equations. + +Conservative variable vector: W = [rho, rho*u, rho*E] + where E = e + u^2/2 is specific total energy, + e = P / (rho * (gamma - 1)) is specific internal energy. + +Physical flux: F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)] + +Available solvers: + - hll_flux: HLL (Harten-Lax-van Leer) two-wave approximate solver + - roe_flux: Roe linearized solver with Harten-Hyman entropy fix +""" + +import numpy as np + + +# --------------------------------------------------------------------------- +# Helper: recover primitives + physical flux from a conservative state +# --------------------------------------------------------------------------- +def _primitives(W, gamma): + """Return (rho, u, P, a, H) from conservative W = [rho, rho*u, rho*E].""" + rho = W[0] + if rho <= 0: + raise ValueError(f"Non-positive density: rho={rho}, W={W}") + u = W[1] / rho + E = W[2] + P = (gamma - 1) * (E - 0.5 * rho * u ** 2) + if P <= 0: + raise ValueError(f"Non-positive pressure: P={P}, W={W}") + a = np.sqrt(gamma * P / rho) + H = (E + P) / rho # specific total enthalpy + return rho, u, P, a, H + + +def _physical_flux(rho, u, P, E): + """Physical Euler flux from primitives + total energy density.""" + return np.array([ + rho * u, + rho * u ** 2 + P, + u * (E + P), + ]) + + +# --------------------------------------------------------------------------- +# HLL solver +# --------------------------------------------------------------------------- +def hll_flux(W_L, W_R, gamma): + """ + Compute the HLL numerical flux at the interface between two states. + + Parameters + ---------- + W_L, W_R : array-like of shape (3,) + Left and right conservative state vectors. + gamma : float + Ratio of specific heats. + + Returns + ------- + np.ndarray of shape (3,) + HLL numerical flux vector. + + Raises + ------ + ValueError + If either state has non-positive density or pressure. + """ + rho_L, u_L, p_L, a_L, _ = _primitives(W_L, gamma) + rho_R, u_R, p_R, a_R, _ = _primitives(W_R, gamma) + E_L, E_R = W_L[2], W_R[2] + + F_L = _physical_flux(rho_L, u_L, p_L, E_L) + F_R = _physical_flux(rho_R, u_R, p_R, E_R) + + # --- HLL wave-speed estimates (Davis) --- + S_L = min(u_L - a_L, u_R - a_R) + S_R = max(u_L + a_L, u_R + a_R) + + # --- HLL flux, piecewise on wave configuration --- + if S_L >= 0: + return F_L + if S_R <= 0: + return F_R + W_L_arr = np.asarray(W_L, dtype=float) + W_R_arr = np.asarray(W_R, dtype=float) + return (S_R * F_L - S_L * F_R + S_L * S_R * (W_R_arr - W_L_arr)) / (S_R - S_L) + + +# --------------------------------------------------------------------------- +# Roe solver with Harten-Hyman entropy fix +# --------------------------------------------------------------------------- +def roe_flux(W_L, W_R, gamma): + """ + Compute the Roe linearized numerical flux with Harten-Hyman entropy fix. + + The Roe solver resolves all three waves (left acoustic, contact/entropy, + right acoustic) and is more accurate than HLL at contact discontinuities. + + Parameters + ---------- + W_L, W_R : array-like of shape (3,) + Left and right conservative state vectors. + gamma : float + Ratio of specific heats. + + Returns + ------- + np.ndarray of shape (3,) + Roe numerical flux vector. + + Raises + ------ + ValueError + If either state has non-positive density or pressure. + """ + rho_L, u_L, p_L, a_L, H_L = _primitives(W_L, gamma) + rho_R, u_R, p_R, a_R, H_R = _primitives(W_R, gamma) + + F_L = _physical_flux(rho_L, u_L, p_L, W_L[2]) + F_R = _physical_flux(rho_R, u_R, p_R, W_R[2]) + + # --- Roe-averaged quantities (density-weighted) --- + sqrt_rL = np.sqrt(rho_L) + sqrt_rR = np.sqrt(rho_R) + denom = sqrt_rL + sqrt_rR + + rho_hat = sqrt_rL * sqrt_rR # geometric mean density + u_hat = (sqrt_rL * u_L + sqrt_rR * u_R) / denom + H_hat = (sqrt_rL * H_L + sqrt_rR * H_R) / denom + a_hat_sq = (gamma - 1) * (H_hat - 0.5 * u_hat ** 2) + if a_hat_sq <= 0: + return hll_flux(W_L, W_R, gamma) # fallback + a_hat = np.sqrt(a_hat_sq) + + # --- Eigenvalues of the Roe matrix --- + lam1 = u_hat - a_hat # left acoustic + lam2 = u_hat # entropy / contact + lam3 = u_hat + a_hat # right acoustic + + # --- Wave strengths (jump decomposition onto eigenvectors) --- + dp = p_R - p_L + du = u_R - u_L + drho = rho_R - rho_L + + alpha_1 = (dp - rho_hat * a_hat * du) / (2.0 * a_hat ** 2) + alpha_2 = drho - dp / (a_hat ** 2) + alpha_3 = (dp + rho_hat * a_hat * du) / (2.0 * a_hat ** 2) + + # --- Right eigenvectors --- + r1 = np.array([1.0, u_hat - a_hat, H_hat - u_hat * a_hat]) + r2 = np.array([1.0, u_hat, 0.5 * u_hat ** 2]) + r3 = np.array([1.0, u_hat + a_hat, H_hat + u_hat * a_hat]) + + # --- Harten-Hyman entropy fix --- + # Prevents unphysical expansion shocks at sonic points + eps1 = max(0.0, lam1 - (u_L - a_L), (u_R - a_R) - lam1) + eps3 = max(0.0, lam3 - (u_L + a_L), (u_R + a_R) - lam3) + + abs_lam1 = abs(lam1) + abs_lam2 = abs(lam2) + abs_lam3 = abs(lam3) + + if abs_lam1 < eps1: + abs_lam1 = (lam1 ** 2 + eps1 ** 2) / (2.0 * eps1) + if abs_lam3 < eps3: + abs_lam3 = (lam3 ** 2 + eps3 ** 2) / (2.0 * eps3) + + # --- Roe flux: F = 0.5*(F_L + F_R) - 0.5 * sum(alpha_k |lam_k| r_k) --- + return 0.5 * (F_L + F_R) - 0.5 * ( + alpha_1 * abs_lam1 * r1 + + alpha_2 * abs_lam2 * r2 + + alpha_3 * abs_lam3 * r3 + ) + + +# --------------------------------------------------------------------------- +# Dispatcher +# --------------------------------------------------------------------------- +_SOLVERS = { + 'hll': hll_flux, + 'roe': roe_flux, +} + + +def get_riemann_solver(name): + """ + Return the Riemann flux function for the given solver name. + + Parameters + ---------- + name : str + Solver name: 'hll' or 'roe'. + + Returns + ------- + callable + A function with signature (W_L, W_R, gamma) -> np.ndarray(3,). + """ + key = name.lower() + if key not in _SOLVERS: + raise ValueError( + f"Unknown Riemann solver '{name}'. Available: {list(_SOLVERS.keys())}" + ) + return _SOLVERS[key] diff --git a/src/solver.py b/src/solver.py new file mode 100644 index 0000000..dba0adf --- /dev/null +++ b/src/solver.py @@ -0,0 +1,115 @@ +""" +Time-loop driver for the 0D-1D coupled tank-pipe simulation. + +Per time step (per spec §3.1): + 1. Compute CFL-limited dt from pipe's max wave speed + 2. Freeze ghost states from current tank states + 3. Compute two boundary HLL fluxes (left and right) + 4. Advance pipe by one step using those two fluxes (pipe.step handles + the internal fluxes itself) + 5. Advance both tanks using the SAME two boundary fluxes * area + -> this "flux doubling" is the mechanism that makes system mass + and energy strictly conserved to machine precision + 6. Advance time + 7. Append snapshot to history +""" +import numpy as np + + +def run(tank1, tank2, pipe, t_end, cfl, verbose=False, log_every=100): + """ + Run the coupled tank-pipe simulation from t=0 to t=t_end. + + Parameters + ---------- + tank1, tank2 : Tank + Upstream and downstream tanks. tank1 connects to pipe.W[:, 0], + tank2 connects to pipe.W[:, -1]. + pipe : Pipe + 1D pipe instance with initial state already set. + t_end : float + End time in seconds. + cfl : float + CFL number in (0, 1]. + verbose : bool, default False + If True, print step-progress info every `log_every` steps. + log_every : int, default 100 + Logging interval when verbose=True. + + Returns + ------- + dict + History with keys 't', 'P1', 'T1', 'P2', 'T2' (all 1D arrays + of shape (n_steps,)), and 'W_hist' of shape (n_steps, 3, N). + """ + history = { + 't': [], + 'P1': [], 'T1': [], + 'P2': [], 'T2': [], + 'W_hist': [], + } + + t = 0.0 + step = 0 + + while t < t_end: + # --- Phase 1: CFL time step --- + a_max = pipe.max_wave_speed() + dt = cfl * pipe.dx / a_max + dt = min(dt, t_end - t) + if dt < 1e-12: + raise RuntimeError( + f"dt degenerate at step {step}: dt={dt:.3e}, a_max={a_max:.3e}" + ) + + # --- Phase 2: freeze tank ghost states (snapshot for this step) --- + W_ghost_L = tank1.ghost_state() + W_ghost_R = tank2.ghost_state() + + # --- Phase 3: two boundary fluxes (solver-level, same solver as pipe) --- + flux_L = pipe._flux_fn(W_ghost_L, pipe.W[:, 0], pipe.gamma) + flux_R = pipe._flux_fn(pipe.W[:, -1], W_ghost_R, pipe.gamma) + + # --- Phase 4: advance pipe (internal fluxes handled inside) --- + pipe.step(flux_L, flux_R, dt) + + # --- Phase 5: advance tanks with the SAME boundary fluxes * area --- + fL_A = flux_L * pipe.area + fR_A = flux_R * pipe.area + # Left boundary flux is "rightward positive"; tank1 loses that mass + tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=dt, sign=-1) + # Right boundary flux is "rightward positive"; tank2 gains that mass + tank2.apply_flux(mdot=fR_A[0], edot=fR_A[2], dt=dt, sign=+1) + + # --- Phase 6: advance time --- + t += dt + step += 1 + + # --- Phase 7: record history --- + history['t'].append(t) + history['P1'].append(tank1.P) + history['T1'].append(tank1.T) + history['P2'].append(tank2.P) + history['T2'].append(tank2.T) + history['W_hist'].append(pipe.W.copy()) + + if verbose and step % log_every == 0: + _, u, _, _ = pipe.primitives() + print( + f"step={step:6d} t={t:.5f} dt={dt:.2e} " + f"P1={tank1.P/1e6:7.4f}MPa P2={tank2.P/1e6:7.4f}MPa " + f"max|u|={float(np.max(np.abs(u))):7.1f}m/s" + ) + + if step == 0: + raise RuntimeError("solver.run() exited without taking any step") + + # Convert lists to arrays for downstream consumers + history['t'] = np.asarray(history['t']) + history['P1'] = np.asarray(history['P1']) + history['T1'] = np.asarray(history['T1']) + history['P2'] = np.asarray(history['P2']) + history['T2'] = np.asarray(history['T2']) + history['W_hist'] = np.stack(history['W_hist']) # shape (n_steps, 3, N) + + return history diff --git a/src/tank.py b/src/tank.py new file mode 100644 index 0000000..ec1dd0f --- /dev/null +++ b/src/tank.py @@ -0,0 +1,82 @@ +# src/tank.py +""" +0D lumped-parameter tank for ideal gas. The tank's *primary* state is +(mass, U) where U is total internal energy in joules. Pressure, temperature, +and density are derived properties computed from (mass, U) on demand, so +they are always consistent with the conservation-law updates. + +Conservation laws (u=0 inside tank): + dm/dt = mdot_in (mass) + dU/dt = Hdot_in = mdot_in * h_t,in (energy, open-system first law) + +where h_t is specific total enthalpy. When the tank couples to a 1D pipe +through the HLL boundary flux, flux[0]*A = mdot and flux[2]*A = Hdot +automatically — see solver.py. +""" + +import numpy as np + + +class Tank: + def __init__(self, V, P_init, T_init, gamma, R_gas): + self.V = V + self.gamma = gamma + self.R = R_gas + rho = P_init / (R_gas * T_init) + self.mass = rho * V + # For u=0, total internal energy equals rho*e*V = P*V / (gamma-1) + self.U = P_init * V / (gamma - 1) + + @property + def rho(self): + return self.mass / self.V + + @property + def T(self): + return (self.U / self.mass) * (self.gamma - 1) / self.R + + @property + def P(self): + return self.rho * self.R * self.T + + def ghost_state(self): + """ + Return the conservative variable vector [rho, rho*u, rho*E] that + represents this tank as a ghost cell for the 1D pipe solver. + Since u_ghost = 0, rho*u = 0 and rho*E = P/(gamma-1). + """ + return np.array([self.rho, 0.0, self.P / (self.gamma - 1)]) + + def apply_flux(self, mdot, edot, dt, sign): + """ + Update (mass, U) from one time step of boundary flux. + + Parameters + ---------- + mdot : float + Mass flux across the interface in kg/s (already multiplied + by pipe cross-sectional area). Sign is the "outward normal" + convention of the pipe: positive = pipe-rightward. + edot : float + Total enthalpy rate in W (= flux[2] * A), same convention. + dt : float + Time-step size in seconds. + sign : int (+1 or -1) + Orientation for this tank. For an upstream tank whose gas + flows "out to the right" into the pipe, the HLL left-boundary + flux has mdot > 0, so sign = -1 (tank loses mass). + For a downstream tank receiving gas from the right boundary + with mdot > 0 entering, sign = +1. + + Raises + ------ + RuntimeError + If the tank's mass becomes non-positive after the update. + """ + self.mass += sign * mdot * dt + self.U += sign * edot * dt + if self.mass <= 0: + raise RuntimeError( + f"Tank mass non-positive after apply_flux: mass={self.mass}, " + f"mdot={mdot}, edot={edot}, dt={dt}, sign={sign}" + ) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/cryo_tank/__init__.py b/tests/cryo_tank/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/cryo_tank/conftest.py b/tests/cryo_tank/conftest.py new file mode 100644 index 0000000..e25beea --- /dev/null +++ b/tests/cryo_tank/conftest.py @@ -0,0 +1,6 @@ +import os +import sys + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_ROOT = os.path.dirname(os.path.dirname(_HERE)) +sys.path.insert(0, os.path.join(_ROOT, "src")) diff --git a/tests/cryo_tank/test_heat_leak.py b/tests/cryo_tank/test_heat_leak.py new file mode 100644 index 0000000..c0fcf28 --- /dev/null +++ b/tests/cryo_tank/test_heat_leak.py @@ -0,0 +1,47 @@ +# tests/cryo_tank/test_heat_leak.py +"""Tests for heat leak models.""" +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 diff --git a/tests/cryo_tank/test_integration.py b/tests/cryo_tank/test_integration.py new file mode 100644 index 0000000..adb720d --- /dev/null +++ b/tests/cryo_tank/test_integration.py @@ -0,0 +1,79 @@ +# tests/cryo_tank/test_integration.py +"""Integration tests for the cryogenic tank simulation.""" +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 diff --git a/tests/cryo_tank/test_properties.py b/tests/cryo_tank/test_properties.py new file mode 100644 index 0000000..984c4f6 --- /dev/null +++ b/tests/cryo_tank/test_properties.py @@ -0,0 +1,67 @@ +# 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) + # 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 diff --git a/tests/cryo_tank/test_tank_model.py b/tests/cryo_tank/test_tank_model.py new file mode 100644 index 0000000..3a821fa --- /dev/null +++ b/tests/cryo_tank/test_tank_model.py @@ -0,0 +1,78 @@ +# tests/cryo_tank/test_tank_model.py +"""Tests for CryoTank model.""" +import pytest +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.305) < 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 diff --git a/tests/test_integration.py b/tests/test_integration.py new file mode 100644 index 0000000..55907ab --- /dev/null +++ b/tests/test_integration.py @@ -0,0 +1,55 @@ +""" +End-to-end conservation tests: after a short run, total system mass +and total system energy must be conserved to machine precision. These +are the load-bearing tests for the "flux doubling" coupling mechanism. +""" +import numpy as np +import pytest +from tank import Tank +from pipe import Pipe +from solver import run + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def _build_scenario(): + """Default blowdown scenario at smaller scale — same physics, short run.""" + tank1 = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + tank2 = Tank(V=10.0, P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=101325.0, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + return tank1, tank2, pipe + + +def _total_mass(tank1, tank2, pipe): + return tank1.mass + tank2.mass + float(np.sum(pipe.W[0, :] * pipe.area * pipe.dx)) + + +def _total_energy(tank1, tank2, pipe): + return tank1.U + tank2.U + float(np.sum(pipe.W[2, :] * pipe.area * pipe.dx)) + + +def test_total_mass_conserved_over_short_run(): + tank1, tank2, pipe = _build_scenario() + m_init = _total_mass(tank1, tank2, pipe) + run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5) + m_final = _total_mass(tank1, tank2, pipe) + rel_err = abs(m_final - m_init) / m_init + assert rel_err < 1e-10, ( + f"Total mass not conserved: m_init={m_init:.6e}, " + f"m_final={m_final:.6e}, rel_err={rel_err:.2e}" + ) + + +def test_total_energy_conserved_over_short_run(): + tank1, tank2, pipe = _build_scenario() + U_init = _total_energy(tank1, tank2, pipe) + run(tank1, tank2, pipe, t_end=1e-3, cfl=0.5) + U_final = _total_energy(tank1, tank2, pipe) + rel_err = abs(U_final - U_init) / U_init + assert rel_err < 1e-10, ( + f"Total energy not conserved: U_init={U_init:.6e}, " + f"U_final={U_final:.6e}, rel_err={rel_err:.2e}" + ) diff --git a/tests/test_pipe.py b/tests/test_pipe.py new file mode 100644 index 0000000..05dfea8 --- /dev/null +++ b/tests/test_pipe.py @@ -0,0 +1,59 @@ +# tests/test_pipe.py +import numpy as np +import pytest +from pipe import Pipe + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def test_pipe_uniform_initialization_all_cells_identical(): + """ + Initial state: uniform P, T, u=0 across all N cells. + All cells should have identical W. primitives() should round-trip + back to u=0 and P=P_init. + """ + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=1e5, T_init=300.0, + gamma=GAMMA, R_gas=R_GAS) + + # All cells identical + for i in range(1, pipe.N): + assert np.allclose(pipe.W[:, i], pipe.W[:, 0], rtol=1e-14) + + rho, u, P, a = pipe.primitives() + rho_expected = 1e5 / (R_GAS * 300.0) + assert np.allclose(u, 0.0) + assert np.allclose(P, 1e5, rtol=1e-10) + assert np.allclose(rho, rho_expected, rtol=1e-10) + + # Sound speed a = sqrt(gamma * P / rho) + a_expected = np.sqrt(GAMMA * 1e5 / rho_expected) + assert np.allclose(a, a_expected, rtol=1e-10) + + # Geometric sanity + assert pipe.dx == 1.0 / 20 + assert pipe.x_centers.shape == (20,) + assert abs(pipe.x_centers[0] - 0.025) < 1e-14 + assert abs(pipe.x_centers[-1] - 0.975) < 1e-14 + + +def test_pipe_uniform_state_plus_matching_boundary_flux_is_static(): + """ + If all cells have identical W (so all internal HLL fluxes are + identical to the physical flux F(W) = [0, P, 0] for u=0 state), + and we pass in boundary fluxes equal to [0, P, 0] as well, then + every difference (flux[i+1] - flux[i]) is zero, so W must not + change after one step. Verify to machine precision. + """ + pipe = Pipe(L=1.0, D=5e-3, N=20, + P_init=1e5, T_init=300.0, + gamma=GAMMA, R_gas=R_GAS) + + W_before = pipe.W.copy() + # Boundary flux matching the static interior: [rho*u, rho*u^2+P, u*(E+P)] + # with u=0 -> [0, P, 0] + flux_boundary = np.array([0.0, 1e5, 0.0]) + pipe.step(flux_boundary, flux_boundary, dt=1e-5) + assert np.allclose(pipe.W, W_before, atol=1e-6, rtol=1e-12) diff --git a/tests/test_riemann.py b/tests/test_riemann.py new file mode 100644 index 0000000..1174469 --- /dev/null +++ b/tests/test_riemann.py @@ -0,0 +1,74 @@ +# tests/test_riemann.py +import numpy as np +from riemann import hll_flux + + +GAMMA = 1.4 + + +def _to_conservative(rho, u, P, gamma): + """(rho, u, P) -> [rho, rho*u, rho*E] where E = e + u^2/2.""" + return np.array([ + rho, + rho * u, + P / (gamma - 1) + 0.5 * rho * u ** 2, + ]) + + +def _physical_flux(W, gamma): + """F(W) = [rho*u, rho*u^2 + P, u*(rho*E + P)].""" + rho = W[0] + u = W[1] / rho + E = W[2] + P = (gamma - 1) * (E - 0.5 * rho * u ** 2) + return np.array([rho * u, rho * u ** 2 + P, u * (E + P)]) + + +def test_hll_identical_states_returns_physical_flux(): + """ + When W_L == W_R, HLL must return the exact physical flux F(W) + with zero numerical dissipation (the (W_R - W_L) term vanishes). + """ + W = _to_conservative(rho=1.2, u=50.0, P=2.5e5, gamma=GAMMA) + F = hll_flux(W, W, GAMMA) + expected = _physical_flux(W, GAMMA) + assert np.allclose(F, expected, rtol=1e-12), f"F={F}, expected={expected}" + + +def test_hll_equal_pressure_equal_energy_gives_exact_pressure_flux(): + """ + Two stationary states (u=0) with same P but different rho: + - Both have the same energy density E = P/(gamma-1), so the HLL + (W_R - W_L)[2] term vanishes -> exact zero energy flux. + - F_L[1] = F_R[1] = P, and W_L[1] = W_R[1] = 0, so the momentum + flux is exactly P. + - The mass flux is NOT exactly zero for HLL (the density jump + triggers the (W_R - W_L)[0] dissipation term) — this is a known + HLL limitation for stationary contact discontinuities. We do + not assert on F[0] here. + """ + W_L = _to_conservative(rho=10.0, u=0.0, P=1e5, gamma=GAMMA) + W_R = _to_conservative(rho=1.0, u=0.0, P=1e5, gamma=GAMMA) + F = hll_flux(W_L, W_R, GAMMA) + assert abs(F[1] - 1e5) < 1e-6, f"momentum flux should equal P=1e5, got {F[1]}" + assert abs(F[2]) < 1e-8, f"energy flux should be exactly 0, got {F[2]}" + + +def test_hll_sod_shock_tube_directional_fluxes_all_positive(): + """ + Classical Sod initial values: + left: (rho, u, P) = (1.0, 0, 1.0) + right: (rho, u, P) = (0.125, 0, 0.1) + The pressure gradient drives flow from left to right, so the + HLL flux at the interface should have all three components + strictly positive: + F[0] > 0 : mass flux rightward + F[1] > 0 : momentum flux rightward + F[2] > 0 : energy flux rightward + """ + W_L = _to_conservative(rho=1.0, u=0.0, P=1.0, gamma=GAMMA) + W_R = _to_conservative(rho=0.125, u=0.0, P=0.1, gamma=GAMMA) + F = hll_flux(W_L, W_R, GAMMA) + assert F[0] > 0, f"expected positive mass flux, got {F[0]}" + assert F[1] > 0, f"expected positive momentum flux, got {F[1]}" + assert F[2] > 0, f"expected positive energy flux, got {F[2]}" diff --git a/tests/test_tank.py b/tests/test_tank.py new file mode 100644 index 0000000..a340723 --- /dev/null +++ b/tests/test_tank.py @@ -0,0 +1,57 @@ +# tests/test_tank.py +import numpy as np +import pytest +from tank import Tank + + +GAMMA = 1.4 +R_GAS = 287.0 + + +def test_tank_initial_state_matches_ideal_gas(): + """ + Given P, T, V construct a Tank; mass and U should match ideal gas: + rho = P / (R * T) + mass = rho * V + U = P * V / (gamma - 1) (since u=0 inside tank) + T = (U / mass) * (gamma - 1) / R (round-trip) + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + rho_expected = 10e6 / (R_GAS * 300.0) + assert abs(tank.rho - rho_expected) < 1e-9 + assert abs(tank.mass - rho_expected * 5.0) < 1e-6 + assert abs(tank.U - 10e6 * 5.0 / (GAMMA - 1)) < 1e-3 + assert abs(tank.P - 10e6) < 1e-3 + assert abs(tank.T - 300.0) < 1e-9 + + +def test_tank_ghost_state_is_stagnation_conservative_vector(): + """ + ghost_state() should return [rho, 0, P/(gamma-1)] + (u_ghost = 0 per spec §2.3, so total energy density equals + internal energy density = P/(gamma-1)). + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + g = tank.ghost_state() + assert g.shape == (3,) + assert abs(g[0] - tank.rho) < 1e-12 + assert g[1] == 0.0 + assert abs(g[2] - 10e6 / (GAMMA - 1)) < 1e-3 + + +def test_tank_apply_flux_outflow_reduces_mass_energy_and_pressure(): + """ + apply_flux(mdot, edot, dt, sign=-1) should subtract mdot*dt from mass + and edot*dt from U. Derived P should decrease correspondingly. + """ + tank = Tank(V=5.0, P_init=10e6, T_init=300.0, gamma=GAMMA, R_gas=R_GAS) + mass_before = tank.mass + U_before = tank.U + P_before = tank.P + mdot = 1.0 # kg/s + edot = 5e5 # J/s (enthalpy rate) + dt = 1e-3 + tank.apply_flux(mdot=mdot, edot=edot, dt=dt, sign=-1) + assert abs(tank.mass - (mass_before - mdot * dt)) < 1e-12 + assert abs(tank.U - (U_before - edot * dt)) < 1e-9 + assert tank.P < P_before