代码仓库移植

src/README.md

@@ -0,0 +1,3 @@
+# src
+
+Source code for pipe system simulation models and utilities.

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"

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

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

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()

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)

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

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

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

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

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()

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 = [
+        "<!doctype html>",
+        "<html lang='en'>",
+        "<head>",
+        "<meta charset='utf-8'>",
+        "<title>Pipe system simulation summary</title>",
+        "<style>",
+        "body { font-family: Arial, sans-serif; margin: 24px; line-height: 1.45; }",
+        "h1, h2 { margin-bottom: 0.3em; }",
+        "table { border-collapse: collapse; width: 100%; margin: 12px 0 24px; }",
+        "th, td { border: 1px solid #ccc; padding: 8px 10px; text-align: left; vertical-align: top; }",
+        "th { background: #f5f5f5; width: 28%; }",
+        "img { max-width: 100%; height: auto; border: 1px solid #ddd; margin: 8px 0 24px; }",
+        "code { background: #f5f5f5; padding: 1px 4px; }",
+        "</style>",
+        "</head>",
+        "<body>",
+        "<h1>Pipe system simulation summary</h1>",
+        f"<p>Generated from <code>results/history.npz</code>. Final simulation time: {t[-1]:.6f} s.</p>",
+    ]
+
+    for title, rows in sections:
+        parts.append(f"<h2>{html.escape(title)}</h2>")
+        parts.append("<table>")
+        for key, value in rows:
+            parts.append(
+                f"<tr><th>{html.escape(str(key))}</th><td>{html.escape(str(value))}</td></tr>"
+            )
+        parts.append("</table>")
+
+    parts.extend([
+        "<h2>Figures</h2>",
+        "<p><img src='tank_pressure.png' alt='Tank pressure history'></p>",
+        "<p><img src='tank_temperature.png' alt='Tank temperature history'></p>",
+        "<p><img src='pipe_final_profiles.png' alt='Final pipe profiles'></p>",
+        "<p>Animation: <a href='pipe_animation.gif'>pipe_animation.gif</a></p>",
+        "</body>",
+        "</html>",
+    ])
+
+    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))

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()}"
+            )

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]

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

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}"
+            )