模型验证补充

AGENTS.md

@@ -1,12 +1,13 @@
 # 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.
+Core simulation code lives in `src/`. The base pipe blowdown model is grouped under `src/tank_pipe/`, the high-pressure gas cylinder model is under `src/cylinder/`, and the cryogenic tank variant is isolated under `src/cryo_tank/`. `src/main.py` remains a compatibility entry point for the tank-pipe simulation. Tests live 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/main.py`: run the 0D-1D tank-pipe simulation through the compatibility entry point.
+- `python3 src/tank_pipe/main.py`: run the 0D-1D tank-pipe simulation from its package entry point.
 - `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.

README.md

@@ -4,12 +4,20 @@ Test repository for pipe system simulation work.
 
 ## Structure
 
+- `src/tank_pipe/` 0D-1D tank-pipe blowdown simulation
+- `src/cylinder/` high-pressure gas cylinder model
+- `src/cryo_tank/` cryogenic LN2 tank simulation
+- `src/main.py` compatibility entry point for the tank-pipe simulation
 - `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
+## Common Commands
 
-Initial repository scaffold created.
+```bash
+python3 src/main.py
+python3 src/tank_pipe/main.py
+python3 src/cryo_tank/main.py
+pytest -q
+```

examples/__init__.py

@@ -0,0 +1 @@
+"""Example system simulations."""

examples/cryo_tank_cylinder_system.py

@@ -0,0 +1,389 @@
+"""
+Coupled cryogenic tank and upstream high-pressure helium cylinder example.
+
+The cryogenic tank still enforces a constant ullage pressure P_work. The
+helium inlet temperature is fixed to cryo_tank.config.T_IN_HE. The resulting
+tank-side helium boundary is passed upstream to the cylinder as the imposed
+mass and energy outflow condition.
+"""
+import os
+import sys
+from dataclasses import dataclass
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.integrate import solve_ivp
+
+_REPO_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+_SRC_DIR = os.path.join(_REPO_ROOT, "src")
+if _SRC_DIR not in sys.path:
+    sys.path.insert(0, _SRC_DIR)
+
+from cryo_tank.config import (  # noqa: E402
+    V_TOTAL, H_TANK, P_WORKING, T_INIT, ULLAGE_FRACTION,
+    MDOT_IN_LN2, T_IN_LN2, MDOT_OUT_LN2,
+    H_CONV_SURFACE, T_ENV, A_TOTAL, T_IN_HE,
+    T_END, RTOL, ATOL,
+)
+from cryo_tank import properties as prop  # noqa: E402
+from cryo_tank.heat_leak import MLIHeatLeak  # noqa: E402
+from cryo_tank.tank_model import CryoTank  # noqa: E402
+from cylinder import HighPressureGasCylinder  # noqa: E402
+
+
+OUTPUT_DIR = os.path.join(_REPO_ROOT, "results", "examples")
+
+
+@dataclass(frozen=True)
+class HeliumBoundary:
+    """Helium inlet boundary passed between the tank and cylinder."""
+
+    mdot_to_tank: float
+    edot_to_tank: float
+    source_pressure: float
+    source_temperature: float
+    source_enthalpy: float
+    boundary_enthalpy: float
+    tank_inlet_pressure: float
+    tank_inlet_temperature: float
+
+
+def build_default_tank():
+    """Build the default cryogenic tank used by the example."""
+    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,
+    )
+
+
+def sync_fixed_he_inlet_boundary(tank):
+    """Apply the fixed tank helium inlet boundary from cryo_tank.config."""
+    boundary_h = prop.he_h(T_IN_HE, tank.P_work)
+    tank.T_in_he = T_IN_HE
+    tank.h_in_he = boundary_h
+    return boundary_h, T_IN_HE
+
+
+def _set_cylinder_conserved_state(cylinder, mass, U):
+    cylinder.mass = mass
+    cylinder.U = U
+    cylinder._update_state()
+
+
+def _heat_terms(tank, info):
+    T_liq = info['T_liq']
+    T_ull = info['T_ull']
+    liquid_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(liquid_level)
+    A_total_current = A_wet + A_dry
+    if A_total_current > 0.0:
+        Q_leak_liq = Q_leak * A_wet / A_total_current
+        Q_leak_ull = Q_leak * A_dry / A_total_current
+    else:
+        Q_leak_liq = 0.0
+        Q_leak_ull = 0.0
+    return Q_liq_to_ull, Q_leak, Q_leak_liq, Q_leak_ull
+
+
+def tank_rates_and_boundary(tank, cylinder, y_tank):
+    """Return tank ODE rates and the coupled helium inlet boundary."""
+    boundary_h, tank_inlet_T = sync_fixed_he_inlet_boundary(tank)
+    info = tank.derive(y_tank)
+    Q_liq_to_ull, _, Q_leak_liq, Q_leak_ull = _heat_terms(tank, info)
+
+    dm_liq_dt = tank.dm_liq_dt
+    dU_liq_dt = tank._liquid_energy_rate(info, Q_liq_to_ull, Q_leak_liq)
+    dT_ull_dt, mdot_he = tank._solve_ullage_temperature_rate(
+        info, Q_liq_to_ull, Q_leak_ull, dU_liq_dt
+    )
+
+    boundary = HeliumBoundary(
+        mdot_to_tank=mdot_he,
+        edot_to_tank=mdot_he * boundary_h,
+        source_pressure=cylinder.P,
+        source_temperature=cylinder.T,
+        source_enthalpy=cylinder.h,
+        boundary_enthalpy=boundary_h,
+        tank_inlet_pressure=tank.P_work,
+        tank_inlet_temperature=tank_inlet_T,
+    )
+    return np.array([dm_liq_dt, dU_liq_dt, dT_ull_dt]), boundary
+
+
+def run_system(t_end=T_END, rtol=RTOL, atol=ATOL, max_step=10.0,
+               tank=None, cylinder=None):
+    """Run the coupled tank-cylinder system.
+
+    Returns a history dict. The cylinder state is included as conserved state
+    variables ``m_cylinder`` and ``U_cylinder`` plus derived pressure,
+    temperature, density, and boundary quantities.
+    """
+    cylinder = HighPressureGasCylinder() if cylinder is None else cylinder
+    tank = build_default_tank() if tank is None else tank
+    sync_fixed_he_inlet_boundary(tank)
+
+    y0 = np.array([
+        *tank.initial_state(),
+        cylinder.mass,
+        cylinder.U,
+    ])
+
+    def rhs(t, y):
+        _set_cylinder_conserved_state(cylinder, y[3], y[4])
+        tank_rates, boundary = tank_rates_and_boundary(tank, cylinder, y[:3])
+        cylinder_mass_rate = -boundary.mdot_to_tank
+        cylinder_energy_rate = -boundary.edot_to_tank
+        return np.array([
+            tank_rates[0],
+            tank_rates[1],
+            tank_rates[2],
+            cylinder_mass_rate,
+            cylinder_energy_rate,
+        ])
+
+    def liquid_empty_event(t, y):
+        return y[0]
+
+    liquid_empty_event.terminal = True
+    liquid_empty_event.direction = -1
+
+    def cylinder_pressure_event(t, y):
+        _set_cylinder_conserved_state(cylinder, y[3], y[4])
+        return cylinder.P - tank.P_work
+
+    cylinder_pressure_event.terminal = True
+    cylinder_pressure_event.direction = -1
+
+    sol = solve_ivp(
+        rhs,
+        [0.0, t_end],
+        y0,
+        method='RK45',
+        rtol=rtol,
+        atol=atol,
+        max_step=max_step,
+        events=[liquid_empty_event, cylinder_pressure_event],
+        dense_output=True,
+    )
+    if not sol.success:
+        raise RuntimeError(f"Coupled solve failed: {sol.message}")
+
+    return post_process_history(tank, cylinder, sol.t, sol.y)
+
+
+def post_process_history(tank, cylinder, t, y):
+    """Compute tank, cylinder, and boundary histories from solver output."""
+    n = len(t)
+    history = {
+        't': t,
+        'm_liq': y[0],
+        'U_liq': y[1],
+        'T_ull': y[2],
+        'm_cylinder': y[3],
+        'U_cylinder': y[4],
+        'T_liq': np.zeros(n),
+        'T_tank_liq': np.zeros(n),
+        'T_tank_ullage': np.zeros(n),
+        'fill_fraction': np.zeros(n),
+        'liquid_level': np.zeros(n),
+        'V_ull': np.zeros(n),
+        'm_He': np.zeros(n),
+        'U_ull': np.zeros(n),
+        'P_tank': np.zeros(n),
+        'mdot_He': np.zeros(n),
+        'mdot_tank_inlet': np.zeros(n),
+        'mdot_cylinder_out': np.zeros(n),
+        'edot_He': np.zeros(n),
+        'P_cylinder': np.zeros(n),
+        'T_cylinder': np.zeros(n),
+        'rho_cylinder': np.zeros(n),
+        'h_cylinder': np.zeros(n),
+        'P_he_source': np.zeros(n),
+        'T_he_source': np.zeros(n),
+        'h_he_source': np.zeros(n),
+        'P_he_boundary': np.zeros(n),
+        'T_he_boundary': np.zeros(n),
+        'h_he_boundary': np.zeros(n),
+        'P_he_tank_inlet': np.zeros(n),
+        'T_he_tank_inlet': np.zeros(n),
+        'pressure_margin': np.zeros(n),
+        'Q_liq_to_ull': np.zeros(n),
+        'Q_leak': np.zeros(n),
+        'Q_leak_liq': np.zeros(n),
+        'Q_leak_ull': np.zeros(n),
+    }
+
+    for i in range(n):
+        _set_cylinder_conserved_state(cylinder, y[3, i], y[4, i])
+        tank_rates, boundary = tank_rates_and_boundary(tank, cylinder, y[:3, i])
+        info = tank.derive(y[:3, i])
+        Q_liq_to_ull, Q_leak, Q_leak_liq, Q_leak_ull = _heat_terms(tank, info)
+
+        history['T_liq'][i] = info['T_liq']
+        history['T_tank_liq'][i] = info['T_liq']
+        history['T_tank_ullage'][i] = info['T_ull']
+        history['fill_fraction'][i] = info['fill_fraction']
+        history['liquid_level'][i] = info['liquid_level']
+        history['V_ull'][i] = info['V_ull']
+        history['m_He'][i] = info['m_He']
+        history['U_ull'][i] = info['U_ull']
+        history['P_tank'][i] = info['P_He']
+        history['mdot_He'][i] = boundary.mdot_to_tank
+        history['mdot_tank_inlet'][i] = boundary.mdot_to_tank
+        history['mdot_cylinder_out'][i] = boundary.mdot_to_tank
+        history['edot_He'][i] = boundary.edot_to_tank
+        history['P_cylinder'][i] = cylinder.P
+        history['T_cylinder'][i] = cylinder.T
+        history['rho_cylinder'][i] = cylinder.rho
+        history['h_cylinder'][i] = cylinder.h
+        history['P_he_source'][i] = boundary.source_pressure
+        history['T_he_source'][i] = boundary.source_temperature
+        history['h_he_source'][i] = boundary.source_enthalpy
+        history['P_he_boundary'][i] = boundary.tank_inlet_pressure
+        history['T_he_boundary'][i] = boundary.tank_inlet_temperature
+        history['h_he_boundary'][i] = boundary.boundary_enthalpy
+        history['P_he_tank_inlet'][i] = boundary.tank_inlet_pressure
+        history['T_he_tank_inlet'][i] = boundary.tank_inlet_temperature
+        history['pressure_margin'][i] = boundary.source_pressure - boundary.tank_inlet_pressure
+        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
+
+    return history
+
+
+def _save_figure(fig, path):
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    fig.tight_layout()
+    fig.savefig(path, dpi=140)
+    plt.close(fig)
+
+
+def plot_requested_outputs(history, output_dir):
+    """Plot tank and cylinder pressure, temperature, and mass-flow histories."""
+    os.makedirs(output_dir, exist_ok=True)
+    t = history['t']
+
+    fig, (ax_p, ax_t) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
+    ax_p.plot(t, history['P_tank'] / 1e6)
+    ax_p.set_ylabel('Tank pressure [MPa]')
+    ax_p.grid(True)
+
+    ax_t.plot(t, history['T_tank_liq'], label='Liquid')
+    ax_t.plot(t, history['T_tank_ullage'], label='Ullage')
+    ax_t.plot(t, history['T_he_tank_inlet'], label='He inlet', linestyle='--')
+    ax_t.set_xlabel('Time [s]')
+    ax_t.set_ylabel('Tank temperature [K]')
+    ax_t.grid(True)
+    ax_t.legend()
+    _save_figure(
+        fig, os.path.join(output_dir, 'tank_pressure_temperature.png')
+    )
+
+    fig, ax = plt.subplots(figsize=(10, 4.5))
+    ax.plot(t, history['mdot_tank_inlet'] * 1000.0)
+    ax.set_xlabel('Time [s]')
+    ax.set_ylabel('Tank inlet He mass flow [g/s]')
+    ax.grid(True)
+    _save_figure(fig, os.path.join(output_dir, 'tank_inlet_mass_flow.png'))
+
+    fig, (ax_p, ax_t) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
+    ax_p.plot(t, history['P_cylinder'] / 1e6)
+    ax_p.set_ylabel('Cylinder pressure [MPa]')
+    ax_p.grid(True)
+
+    ax_t.plot(t, history['T_cylinder'])
+    ax_t.set_xlabel('Time [s]')
+    ax_t.set_ylabel('Cylinder temperature [K]')
+    ax_t.grid(True)
+    _save_figure(
+        fig, os.path.join(output_dir, 'cylinder_pressure_temperature.png')
+    )
+
+    fig, ax = plt.subplots(figsize=(10, 4.5))
+    ax.plot(t, history['mdot_cylinder_out'] * 1000.0)
+    ax.set_xlabel('Time [s]')
+    ax.set_ylabel('Cylinder outlet He mass flow [g/s]')
+    ax.grid(True)
+    _save_figure(fig, os.path.join(output_dir, 'cylinder_mass_flow.png'))
+
+
+def save_history_csv(history, path):
+    """Save a 1D history dictionary as CSV."""
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    names = list(history.keys())
+    data = np.column_stack([np.asarray(history[name]) for name in names])
+    np.savetxt(path, data, delimiter=',', header=','.join(names), comments='')
+
+
+def save_history_npz(history, path):
+    """Save a history dictionary as compressed NPZ."""
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    np.savez_compressed(path, **history)
+
+
+def main():
+    cylinder = HighPressureGasCylinder()
+    tank = build_default_tank()
+
+    print("Coupled cryo tank + upstream He cylinder")
+    print(f"  t_end = {T_END:.3f} s")
+    print(f"  tank P_work = {tank.P_work / 1e6:.4f} MPa")
+    print(f"  cylinder: P = {cylinder.P / 1e6:.4f} MPa, T = {cylinder.T:.2f} K, "
+          f"m = {cylinder.mass:.4f} kg")
+    print()
+
+    history = run_system(tank=tank, cylinder=cylinder)
+
+    csv_path = os.path.join(OUTPUT_DIR, "cryo_tank_cylinder_system.csv")
+    npz_path = os.path.join(OUTPUT_DIR, "cryo_tank_cylinder_system.npz")
+    save_history_csv(history, csv_path)
+    save_history_npz(history, npz_path)
+    plot_requested_outputs(history, OUTPUT_DIR)
+
+    print("Simulation complete:")
+    print(f"  t_final = {history['t'][-1]:.1f} s")
+    print(f"  tank fill_fraction: {history['fill_fraction'][0]:.1%} -> "
+          f"{history['fill_fraction'][-1]:.1%}")
+    print(f"  tank P: {history['P_tank'][0] / 1e6:.4f} -> "
+          f"{history['P_tank'][-1] / 1e6:.4f} MPa")
+    print(f"  tank T_liq: {history['T_tank_liq'][0]:.2f} -> "
+          f"{history['T_tank_liq'][-1]:.2f} K")
+    print(f"  tank T_ull: {history['T_tank_ullage'][0]:.2f} -> "
+          f"{history['T_tank_ullage'][-1]:.2f} K")
+    print(f"  tank inlet He mdot: {history['mdot_tank_inlet'][0] * 1000:.5f} -> "
+          f"{history['mdot_tank_inlet'][-1] * 1000:.5f} g/s")
+    print(f"  cylinder P: {history['P_cylinder'][0] / 1e6:.4f} -> "
+          f"{history['P_cylinder'][-1] / 1e6:.4f} MPa")
+    print(f"  cylinder T: {history['T_cylinder'][0]:.2f} -> "
+          f"{history['T_cylinder'][-1]:.2f} K")
+    print(f"  cylinder outlet He mdot: {history['mdot_cylinder_out'][0] * 1000:.5f} -> "
+          f"{history['mdot_cylinder_out'][-1] * 1000:.5f} g/s")
+    print(f"  cylinder mass: {history['m_cylinder'][0]:.4f} -> "
+          f"{history['m_cylinder'][-1]:.4f} kg")
+    print(f"  tank-side He inlet T: {history['T_he_tank_inlet'][0]:.2f} -> "
+          f"{history['T_he_tank_inlet'][-1]:.2f} K")
+    print(f"  fixed He boundary h: {history['h_he_boundary'][0]:.2f} -> "
+          f"{history['h_he_boundary'][-1]:.2f} J/kg")
+    print(f"\nOutputs written to {OUTPUT_DIR}/")
+
+
+if __name__ == "__main__":
+    main()

src/README.md

@@ -1,3 +1,20 @@
 # src
 
 Source code for pipe system simulation models and utilities.
+
+## Layout
+
+- `main.py`: compatibility entry point for the tank-pipe simulation.
+- `tank_pipe/`: 0D-1D tank-pipe blowdown model, including configuration, tank, pipe, Riemann solver, time solver, and output helpers.
+- `cylinder/`: high-pressure gas cylinder model.
+- `cryo_tank/`: cryogenic LN2 tank model with helium pressurization, heat leak, solver, output helpers, and module README.
+
+## Entry Points
+
+Run from the repository root:
+
+```bash
+python3 src/main.py
+python3 src/tank_pipe/main.py
+python3 src/cryo_tank/main.py
+```

src/cryo_tank/Index.md

@@ -0,0 +1,710 @@
+## 1. 模型假设
+
+当前模型是液氮贮箱的两区集总参数模型：
+
+- 液相区为液氮。
+- 气枕区为纯氦气。
+- 气枕压力固定为工作压力 `P_work`。
+- 不考虑氮气蒸发、冷凝和气液相间质量传递。
+- 氦气质量不是独立状态量，而是维持恒压所需的派生量。
+- 氦气只允许补入，不允许倒流；若计算得到负的氦气流量，则报告流量钳制为 0。
+
+对应代码：
+
+- 状态和派生量：`src/cryo_tank/tank_model.py:98-130`
+- ODE 右端项：`src/cryo_tank/tank_model.py:132-162`
+- 液相温度、气枕体积和氦气补气方程：`src/cryo_tank/tank_model.py:172-257`
+
+## 2. 状态向量
+
+ODE 状态向量为：
+
+$$
+\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right]
+$$
+
+其中：
+
+$$
+\begin{aligned}
+m_{\mathrm{liq}} &: \text{液氮质量}\ [\mathrm{kg}] \\
+U_{\mathrm{liq}} &: \text{液相总内能}\ [\mathrm{J}] \\
+T_{\mathrm{ull}} &: \text{气枕区温度}\ [\mathrm{K}]
+\end{aligned}
+$$
+
+式中，`m_liq` 为液氮质量，`U_liq` 为液相总内能，`T_ull` 为氦气气枕温度。`m_liq` 和 `U_liq` 是液相守恒方程的状态量，`T_ull` 是气枕能量方程的状态量。
+
+代码位置：`src/cryo_tank/tank_model.py:1-7`，`src/cryo_tank/tank_model.py:87-89`
+
+## 3. 初始状态方程
+
+初始液相体积：
+
+$$
+V_{\mathrm{liq},0}=\left(1-f_{\mathrm{ull},0}\right)V_{\mathrm{total}}
+$$
+
+式中，`V_liq_0` 为初始液相体积，`ullage_fraction` 为初始气枕体积分数，`V_total` 为贮箱总容积。
+
+初始气枕体积：
+
+$$
+V_{\mathrm{ull},0}=f_{\mathrm{ull},0}V_{\mathrm{total}}
+$$
+
+式中，`V_ull_0` 为初始气枕体积。
+
+初始液氮质量：
+
+$$
+\begin{aligned}
+\rho_{\mathrm{liq},0} &= \rho_{\mathrm{LN2}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right) \\
+m_{\mathrm{liq},0} &= \rho_{\mathrm{liq},0}V_{\mathrm{liq},0}
+\end{aligned}
+$$
+
+式中，`rho_liq_0` 为初始液氮密度，`rho_LN2(T, P)` 为液氮在温度 `T`、压力 `P` 下的密度物性函数，`T_init` 为初始温度，`P_work` 为工作压力，`m_liq_0` 为初始液氮质量。
+
+初始液相总内能：
+
+$$
+U_{\mathrm{liq},0}=m_{\mathrm{liq},0}u_{\mathrm{LN2}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right)
+$$
+
+式中，`U_liq_0` 为初始液相总内能，`u_LN2(T, P)` 为液氮比内能。
+
+初始气枕温度：
+
+$$
+T_{\mathrm{ull},0}=T_{\mathrm{init}}
+$$
+
+式中，`T_ull_0` 为初始气枕温度。
+
+初始氦气质量是派生量：
+
+$$
+m_{\mathrm{He},0}=\rho_{\mathrm{He}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right)V_{\mathrm{ull},0}
+$$
+
+式中，`m_He_0` 为初始氦气质量，`rho_He(T, P)` 为氦气密度物性函数。
+
+代码位置：`src/cryo_tank/tank_model.py:71-85`
+
+## 4. 几何方程
+
+横截面积：
+
+$$
+A_{\mathrm{cross}}=\frac{V_{\mathrm{total}}}{H_{\mathrm{tank}}}
+$$
+
+式中，`A_cross` 为圆柱贮箱横截面积，`H_tank` 为贮箱高度。
+
+圆柱直径：
+
+$$
+D=\sqrt{\frac{4A_{\mathrm{cross}}}{\pi}}
+$$
+
+式中，`D` 为等效圆柱直径，`pi` 为圆周率。
+
+侧面积：
+
+$$
+A_{\mathrm{side}}=\pi D H_{\mathrm{tank}}
+$$
+
+式中，`A_side` 为圆柱侧壁面积。
+
+端盖面积：
+
+$$
+A_{\mathrm{cap}}=A_{\mathrm{cross}}
+$$
+
+式中，`A_cap` 为单个端盖面积。
+
+总面积：
+
+$$
+A_{\mathrm{total}}=A_{\mathrm{side}}+2A_{\mathrm{cap}}
+$$
+
+式中，`A_total` 为贮箱外表面总面积，包含侧壁和两个端盖。
+
+代码位置：`src/cryo_tank/tank_model.py:41-48`
+
+液位对应的湿壁和干壁面积：
+
+$$
+\begin{aligned}
+L &= \operatorname{clamp}\left(H_{\mathrm{liq}},0,H_{\mathrm{tank}}\right) \\
+A_{\mathrm{wet}} &= A_{\mathrm{cap}}+\pi D L \\
+A_{\mathrm{dry}} &= A_{\mathrm{cap}}+\pi D\left(H_{\mathrm{tank}}-L\right)
+\end{aligned}
+$$
+
+式中，`liquid_level` 为由液相体积换算得到的液位，`level` 为限制在 `[0, H_tank]` 范围内的液位，`A_wet` 为与液相接触的湿壁面积，`A_dry` 为与气枕接触的干壁面积。
+
+代码位置：`src/cryo_tank/tank_model.py:91-96`
+
+## 5. 派生量方程
+
+给定状态：
+
+$$
+\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right]
+$$
+
+### 5.1 液相派生量
+
+$$
+\begin{aligned}
+u_{\mathrm{liq}} &= \frac{U_{\mathrm{liq}}}{m_{\mathrm{liq}}} \\
+T_{\mathrm{liq}} &= T_{\mathrm{LN2}}\left(u_{\mathrm{liq}},P_{\mathrm{work}}\right) \\
+\rho_{\mathrm{liq}} &= \rho_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right) \\
+V_{\mathrm{liq}} &= \frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}} \\
+H_{\mathrm{liq}} &= \frac{V_{\mathrm{liq}}}{A_{\mathrm{cross}}} \\
+f_{\mathrm{fill}} &= \frac{H_{\mathrm{liq}}}{H_{\mathrm{tank}}}
+\end{aligned}
+$$
+
+式中，`u_liq` 为液氮比内能，`T_liq` 为由比内能和压力反算得到的液相温度，`rho_liq` 为当前液氮密度，`V_liq` 为液相体积，`liquid_level` 为液位高度，`fill_fraction` 为液位高度占贮箱高度的比例。
+
+代码位置：`src/cryo_tank/tank_model.py:107-113`
+
+### 5.2 气枕派生量
+
+$$
+\begin{aligned}
+V_{\mathrm{ull}} &= V_{\mathrm{total}}-V_{\mathrm{liq}} \\
+P_{\mathrm{He}} &= P_{\mathrm{work}} \\
+\rho_{\mathrm{He}} &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{He}}\right) \\
+m_{\mathrm{He}} &= \rho_{\mathrm{He}}V_{\mathrm{ull}} \\
+U_{\mathrm{ull}} &= m_{\mathrm{He}}u_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{He}}\right)
+\end{aligned}
+$$
+
+式中，`V_ull` 为气枕体积，`P_He` 为氦气压力，`rho_He` 为当前氦气密度，`m_He` 为维持恒压所需的氦气质量，`U_ull` 为气枕氦气总内能，`u_He(T, P)` 为氦气比内能。
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+## 6. 换热方程
+
+### 6.1 液相与气枕界面换热
+
+$$
+Q_{\mathrm{liq}\to\mathrm{ull}}=h_{\mathrm{conv}}A_{\mathrm{cross}}\left(T_{\mathrm{liq}}-T_{\mathrm{ull}}\right)
+$$
+
+式中，`Q_liq_to_ull` 为液相到气枕的界面换热率，`h_conv` 为界面对流换热系数，`A_cross` 为气液界面面积，`T_liq` 和 `T_ull` 分别为液相温度和气枕温度。
+
+符号约定：
+
+$$
+\begin{aligned}
+Q_{\mathrm{liq}\to\mathrm{ull}} &> 0 && \text{热量从液相传给气枕} \\
+Q_{\mathrm{liq}\to\mathrm{ull}} &< 0 && \text{热量从气枕传给液相}
+\end{aligned}
+$$
+
+代码位置：`src/cryo_tank/tank_model.py:142-143`
+
+### 6.2 外界总漏热
+
+当前主程序默认使用 MLI 漏热模型：
+
+$$
+Q_{\mathrm{leak}}=A_{\mathrm{total}}q_{\mathrm{MLI}}
+$$
+
+式中，`Q_leak` 为外界向贮箱的总漏热率，`q_mli` 为 MLI 单位面积漏热热流密度。
+
+代码位置：
+
+- 默认创建：`src/cryo_tank/main.py:32-33`
+- MLI 方程：`src/cryo_tank/heat_leak.py:21-40`
+
+代码中还提供 Foam 漏热模型：
+
+$$
+\begin{aligned}
+T_{\mathrm{mean}} &= \frac{T_{\mathrm{inner}}+T_{\mathrm{env}}}{2} \\
+Q_{\mathrm{leak}} &= A_{\mathrm{total}}k_{\mathrm{eff}}\left(T_{\mathrm{mean}}\right)\frac{T_{\mathrm{env}}-T_{\mathrm{inner}}}{\delta}
+\end{aligned}
+$$
+
+式中，`T_mean` 为泡沫层平均温度，`T_inner` 为贮箱内侧温度，`T_env` 为环境温度，`k_eff` 为有效导热系数，`delta` 为泡沫保温层厚度。
+
+如果 `k_eff` 是常数，则直接使用该常数。
+
+代码位置：`src/cryo_tank/heat_leak.py:43-70`
+
+### 6.3 漏热分配
+
+总漏热按湿壁面积和干壁面积分配：
+
+$$
+\begin{aligned}
+Q_{\mathrm{leak,liq}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{wet}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} \\
+Q_{\mathrm{leak,ull}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{dry}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}}
+\end{aligned}
+$$
+
+式中，`Q_leak_liq` 为分配到液相的漏热，`Q_leak_ull` 为分配到气枕的漏热。
+
+代码位置：`src/cryo_tank/tank_model.py:145-149`
+
+## 7. 液氮质量方程
+
+液氮质量变化率：
+
+$$
+\frac{dm_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}
+$$
+
+式中，`dm_liq/dt` 为液氮质量变化率，`mdot_in_ln2` 为液氮入口质量流量，`mdot_out_ln2` 为液氮出口质量流量。
+
+代码位置：
+
+- 净流量预计算：`src/cryo_tank/tank_model.py:64-65`
+- ODE 中使用：`src/cryo_tank/tank_model.py:151-152`
+
+## 8. 液相能量方程
+
+入口液氮焓：
+
+$$
+h_{\mathrm{in,LN2}}=h_{\mathrm{LN2}}\left(T_{\mathrm{in,LN2}},P_{\mathrm{work}}\right)
+$$
+
+式中，`h_in_ln2` 为入口液氮比焓，`h_LN2(T, P)` 为液氮比焓物性函数，`T_in_ln2` 为入口液氮温度。
+
+当前液相出口焓：
+
+$$
+h_{\mathrm{liq}}=h_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right)
+$$
+
+式中，`h_liq` 为当前液相出口比焓，模型假定出口液氮状态等于贮箱液相主体状态。
+
+液相总内能变化率：
+
+$$
+\frac{dU_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}h_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}h_{\mathrm{liq}}-Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,liq}}
+$$
+
+式中，`dU_liq/dt` 为液相总内能变化率；入口液氮携带焓流 `mdot_in_ln2 * h_in_ln2`，出口液氮带走焓流 `mdot_out_ln2 * h_liq`，液相向气枕换热为 `Q_liq_to_ull`，外界向液相漏热为 `Q_leak_liq`。
+
+代码位置：
+
+- 入口焓预计算：`src/cryo_tank/tank_model.py:60-62`
+- 能量方程调用：`src/cryo_tank/tank_model.py:153-155`
+- 能量方程实现：`src/cryo_tank/tank_model.py:164-170`
+
+## 9. 气枕温度和氦气质量导数方程
+
+当前模型把 `T_ull` 作为 ODE 状态。氦气质量由恒压关系派生：
+
+$$
+m_{\mathrm{He}}=\rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right)V_{\mathrm{ull}}
+$$
+
+式中，`m_He` 是派生氦气质量，不是 ODE 独立状态。
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+### 9.1 气枕体积变化率
+
+液相体积由液氮质量和恒压液氮密度决定：
+
+$$
+V_{\mathrm{liq}}=\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right)}
+$$
+
+气枕体积为：
+
+$$
+V_{\mathrm{ull}}=V_{\mathrm{total}}-V_{\mathrm{liq}}
+$$
+
+对液相体积求全导数：
+
+$$
+\frac{dV_{\mathrm{liq}}}{dt}=\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}-\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2}
+\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt}
+$$
+
+其中液相温度变化率由 `u_liq = U_liq / m_liq` 得到：
+
+$$
+\frac{dT_{\mathrm{liq}}}{dt}=\frac{\frac{dU_{\mathrm{liq}}}{dt}-u_{\mathrm{liq}}\frac{dm_{\mathrm{liq}}}{dt}}{m_{\mathrm{liq}}\left(\frac{\partial u_{\mathrm{liq}}}{\partial T}\right)_P}
+$$
+
+因此气枕体积变化率为：
+
+$$
+\frac{dV_{\mathrm{ull}}}{dt}=-\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}+\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2}
+\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt}
+$$
+
+式中，`dV_ull/dt` 为气枕体积变化率，第一项为液氮质量变化导致的体积变化，第二项为液氮温度变化导致密度变化后产生的体积修正项。液氮的恒压密度温度导数和恒压内能温度导数由 CoolProp 提供。
+
+代码位置：
+
+- LN2 恒压偏导：`src/cryo_tank/properties.py:47-56`
+- 液相温度变化率：`src/cryo_tank/tank_model.py:172-180`
+- 气枕体积变化率：`src/cryo_tank/tank_model.py:182-194`
+
+### 9.2 气枕区非氦气入口热量
+
+$$
+Q_{\mathrm{ull,noHe}}=Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,ull}}
+$$
+
+式中，`Q_ull_no_he` 为不含氦气入口焓流的气枕净受热率，包括液相传入气枕的界面换热和分配到气枕干壁的外界漏热。
+
+代码位置：`src/cryo_tank/tank_model.py:221`
+
+### 9.3 恒压氦气状态函数
+
+在固定压力 `P_work` 下，给定 `T_ull` 和 `V_ull`：
+
+$$
+\begin{aligned}
+\rho &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right) \\
+m(T,V) &= \rho V \\
+U(T,V) &= m(T,V)u_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right)
+\end{aligned}
+$$
+
+式中，`rho` 为氦气密度，`T` 表示 `T_ull`，`V` 表示 `V_ull`，`m(T, V)` 为恒压下由温度和体积确定的氦气质量，`U(T, V)` 为恒压下由温度和体积确定的气枕总内能。
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+### 9.4 解析偏导数
+
+代码不再对温度和体积做有限差分，而是使用恒压物性导数：
+
+$$
+\begin{aligned}
+\left(\frac{\partial\rho}{\partial T}\right)_P &= \rho_T\big|_P \\
+\left(\frac{\partial u}{\partial T}\right)_P &= u_T\big|_P
+\end{aligned}
+$$
+
+式中，`drhodT_P` 为恒压氦气密度对温度的偏导数，`dudT_P` 为恒压氦气比内能对温度的偏导数。
+
+CoolProp 偏导接口位置：`src/cryo_tank/properties.py:113-122`
+
+由此得到：
+
+$$
+\begin{aligned}
+m_T &= V_{\mathrm{ull}}\left(\frac{\partial\rho}{\partial T}\right)_P \\
+m_V &= \rho \\
+U_T &= V_{\mathrm{ull}}\left[u\left(\frac{\partial\rho}{\partial T}\right)_P+\rho\left(\frac{\partial u}{\partial T}\right)_P\right] \\
+U_V &= \rho u
+\end{aligned}
+$$
+
+式中，`m_T = (partial m / partial T)_V`，`m_V = (partial m / partial V)_T`，`U_T = (partial U / partial T)_V`，`U_V = (partial U / partial V)_T`。这些偏导用于把 `dm_He/dt` 和 `dU_ull/dt` 展开为 `dT_ull/dt` 与 `dV_ull/dt` 的线性组合。
+
+推导过程如下。恒压下 `rho = rho(T, P_work)`，`u = u(T, P_work)`，所以：
+
+$$
+m(T,V)=\rho\left(T,P_{\mathrm{work}}\right)V
+$$
+
+对温度和体积求偏导：
+
+$$
+\begin{aligned}
+m_T &= \left(\frac{\partial m}{\partial T}\right)_V
+     = V\left(\frac{\partial\rho}{\partial T}\right)_P
+     = V_{\mathrm{ull}}\rho_T\big|_P \\
+m_V &= \left(\frac{\partial m}{\partial V}\right)_T
+     = \rho
+\end{aligned}
+$$
+
+气枕总内能为：
+
+$$
+U(T,V)=m(T,V)u\left(T,P_{\mathrm{work}}\right)=\rho\left(T,P_{\mathrm{work}}\right)Vu\left(T,P_{\mathrm{work}}\right)
+$$
+
+对温度和体积求偏导：
+
+$$
+\begin{aligned}
+U_T &= \left(\frac{\partial U}{\partial T}\right)_V
+     = V\left[u\left(\frac{\partial\rho}{\partial T}\right)_P+\rho\left(\frac{\partial u}{\partial T}\right)_P\right] \\
+    &= V_{\mathrm{ull}}\left(u\rho_T\big|_P+\rho u_T\big|_P\right) \\
+U_V &= \left(\frac{\partial U}{\partial V}\right)_T
+     = \rho u
+\end{aligned}
+$$
+
+代码位置：`src/cryo_tank/tank_model.py:196-209`
+
+### 9.5 气枕温度变化率
+
+氦气入口焓：
+
+$$
+h_{\mathrm{in,He}}=h_{\mathrm{He}}\left(T_{\mathrm{in,He}},P_{\mathrm{work}}\right)
+$$
+
+式中，`h_in_he` 为入口氦气比焓，`h_He(T, P)` 为氦气比焓物性函数，`T_in_he` 为入口氦气温度。
+
+氦气气枕能量方程采用开口、变体积控制体的一阶热力学形式：
+
+$$
+\frac{dU_{\mathrm{ull}}}{dt}=\dot{m}_{\mathrm{He}}h_{\mathrm{in,He}}+Q_{\mathrm{ull,noHe}}-P_{\mathrm{work}}\frac{dV_{\mathrm{ull}}}{dt}
+$$
+
+式中，`dU_ull/dt` 为气枕氦气总内能变化率，`mdot_He * h_in_he` 为补入氦气带入的焓流，`Q_ull_no_he` 为界面换热与外界漏热给气枕的热输入，`P_work * dV_ull/dt` 为气枕边界膨胀功项。当前模型无气枕出口流量，因此没有出口焓流项。
+
+由恒压状态函数的全微分：
+
+$$
+\begin{aligned}
+\frac{dm_{\mathrm{He}}}{dt} &= m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt} \\
+\frac{dU_{\mathrm{ull}}}{dt} &= U_T\frac{dT_{\mathrm{ull}}}{dt}+U_V\frac{dV_{\mathrm{ull}}}{dt}
+\end{aligned}
+$$
+
+把质量导数代入气枕能量方程：
+
+$$
+U_T\frac{dT_{\mathrm{ull}}}{dt}+U_V\frac{dV_{\mathrm{ull}}}{dt}=h_{\mathrm{in,He}}\left(m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt}\right)+Q_{\mathrm{ull,noHe}}-P_{\mathrm{work}}\frac{dV_{\mathrm{ull}}}{dt}
+$$
+
+整理 `dT_ull/dt` 项：
+
+$$
+\left(U_T-h_{\mathrm{in,He}}m_T\right)\frac{dT_{\mathrm{ull}}}{dt}=Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt}
+$$
+
+气枕温度变化率：
+
+$$
+\frac{dT_{\mathrm{ull}}}{dt}=\frac{Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt}}{U_T-h_{\mathrm{in,He}}m_T}
+$$
+
+式中，分子表示扣除体积变化、边界功和入口质量变化耦合项后的气枕有效热输入，分母表示在恒压约束和入口补气耦合下气枕对温度变化的等效热容项。
+
+代码位置：`src/cryo_tank/tank_model.py:223-230`
+
+若分母过小：
+
+$$
+\left|U_T-h_{\mathrm{in,He}}m_T\right|<10^{-30}
+\quad\Longrightarrow\quad
+\begin{cases}
+\dfrac{dT_{\mathrm{ull}}}{dt}=0 \\
+\dot{m}_{\mathrm{He}}=0
+\end{cases}
+$$
+
+式中，`1e-30` 是数值保护阈值，用于避免除以过小分母。
+
+代码位置：`src/cryo_tank/tank_model.py:223-226`
+
+### 9.6 氦气质量变化率
+
+$$
+\frac{dm_{\mathrm{He}}}{dt}=\dot{m}_{\mathrm{He}}
+$$
+
+式中，`dm_He/dt` 为恒压约束下氦气质量变化率，`mdot_He` 为需要补入的氦气质量流量。
+
+其中：
+
+$$
+\dot{m}_{\mathrm{He}}=m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt}
+$$
+
+式中，第一项为气枕温度变化导致的氦气质量变化，第二项为气枕体积变化导致的氦气质量变化。
+
+代码位置：`src/cryo_tank/tank_model.py:231`
+
+若 `mdot_He < 0`，报告流量钳制为 0：
+
+$$
+\dot{m}_{\mathrm{He}}=0
+$$
+
+式中，钳制只作用于报告的补气流量；当前 ODE 状态 `T_ull` 的导数仍由未钳制前的恒压能量方程求得。
+
+代码位置：`src/cryo_tank/tank_model.py:233-238`
+
+## 10. 完整 ODE 方程组
+
+状态向量：
+
+$$
+\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right]
+$$
+
+式中，`y` 为 ODE 状态向量，三个分量分别为液氮质量、液相总内能和气枕温度。
+
+ODE 方程组：
+
+$$
+\frac{dm_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}
+$$
+
+式中，液氮质量变化率等于入口液氮质量流量减去出口液氮质量流量。
+
+$$
+\frac{dU_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}h_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}h_{\mathrm{liq}}-Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,liq}}
+$$
+
+式中，液相总内能变化率由入口焓流、出口焓流、液相向气枕换热和液相漏热共同决定。
+
+$$
+\frac{dT_{\mathrm{ull}}}{dt}=\frac{Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt}}{U_T-h_{\mathrm{in,He}}m_T}
+$$
+
+式中，气枕温度变化率由氦气气枕能量方程和恒压状态方程联立得到。
+
+其中：
+
+$$
+\begin{aligned}
+h_{\mathrm{liq}} &= h_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right) \\
+Q_{\mathrm{liq}\to\mathrm{ull}} &= h_{\mathrm{conv}}A_{\mathrm{cross}}\left(T_{\mathrm{liq}}-T_{\mathrm{ull}}\right) \\
+Q_{\mathrm{ull,noHe}} &= Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,ull}} \\
+Q_{\mathrm{leak,liq}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{wet}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} \\
+Q_{\mathrm{leak,ull}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{dry}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}}
+\end{aligned}
+$$
+
+气枕体积变化率使用完整链式法则：
+
+$$
+\begin{aligned}
+\frac{dT_{\mathrm{liq}}}{dt} &= \frac{\frac{dU_{\mathrm{liq}}}{dt}-u_{\mathrm{liq}}\frac{dm_{\mathrm{liq}}}{dt}}{m_{\mathrm{liq}}\left(\frac{\partial u_{\mathrm{liq}}}{\partial T}\right)_P} \\
+\frac{dV_{\mathrm{ull}}}{dt} &= -\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}+\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2}
+\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt}
+\end{aligned}
+$$
+
+式中，`h_liq` 为当前液相比焓，`Q_liq_to_ull` 为液相到气枕界面换热，`Q_ull_no_he` 为不含入口氦气焓流的气枕热输入，`Q_leak_liq` 和 `Q_leak_ull` 分别为分配到液相和气枕的漏热；`dV_ull/dt` 同时包含液氮质量变化和液氮密度随温度变化的贡献。
+
+氦气质量和质量流量是派生量：
+
+$$
+\begin{aligned}
+m_{\mathrm{He}} &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right)V_{\mathrm{ull}} \\
+\dot{m}_{\mathrm{He}} &= m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt}
+\end{aligned}
+$$
+
+式中，`m_He` 为恒压状态方程给出的当前氦气质量，`mdot_He` 为由质量全微分得到的所需补气流量。
+
+代码位置：`src/cryo_tank/tank_model.py:132-257`
+
+## 11. ODE 求解逻辑
+
+求解器入口：
+
+```text
+run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0)
+```
+
+式中，`tank` 为 `CryoTank` 模型实例，`t_end` 为仿真终止时间，`rtol` 和 `atol` 为相对、绝对误差容限，`max_step` 为求解器最大时间步长。
+
+代码位置：`src/cryo_tank/solver.py:22-24`
+
+求解初值：
+
+$$
+\mathbf{y}_0=\operatorname{initial\_state}\left(\mathrm{tank}\right)
+$$
+
+式中，`y0` 为初始 ODE 状态向量。
+
+代码位置：`src/cryo_tank/solver.py:43`
+
+调用 `solve_ivp`：
+
+```text
+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,
+)
+```
+
+式中，`tank.rhs` 为 ODE 右端函数，`[0.0, t_end]` 为积分时间区间，`method="RK45"` 指四、五阶 Runge-Kutta 自适应算法，`events` 用于注册液体排空终止事件，`dense_output=True` 表示生成连续插值解。
+
+代码位置：`src/cryo_tank/solver.py:45-55`
+
+液体排空事件：
+
+$$
+g(t,\mathbf{y})=m_{\mathrm{liq}}
+$$
+
+式中，`event(t, y)` 为事件函数；当状态中的 `m_liq` 到达零时，事件函数到达零。
+
+事件属性：
+
+$$
+\begin{aligned}
+\mathrm{terminal} &= \mathrm{True} \\
+\mathrm{direction} &= -1
+\end{aligned}
+$$
+
+式中，`terminal=True` 表示事件触发后终止积分，`direction=-1` 表示只检测事件函数从正值下降到零的穿越。
+
+含义：
+
+- 当液氮质量下降到 0 时终止求解。
+- 只检测从正到零的方向。
+
+代码位置：`src/cryo_tank/solver.py:14-19`
+
+求解失败时抛出错误；若液体提前排空，则发出警告。
+
+代码位置：`src/cryo_tank/solver.py:57-64`
+
+## 12. 后处理逻辑
+
+求解完成后，`solver.run()` 对每个输出时刻执行：
+
+$$
+\mathrm{info}_i=\operatorname{derive}\left(\mathbf{y}_i\right)
+$$
+
+式中，`y_i` 为某一输出时刻的状态向量，`info` 为根据该状态计算出的温度、体积、质量和换热等派生量字典。
+
+并重新计算：
+
+$$
+\left\{
+Q_{\mathrm{liq}\to\mathrm{ull}},\
+Q_{\mathrm{leak}},\
+Q_{\mathrm{leak,liq}},\
+Q_{\mathrm{leak,ull}},\
+\dot{m}_{\mathrm{He}}
+\right\}
+$$
+
+式中，这些量分别为界面换热、总漏热、液相漏热、气枕漏热和派生氦气补气质量流量，用于输出和结果分析。
+
+其中 `m_He`、`U_ull`、`mdot_He` 均为后处理派生结果，不是独立 ODE 状态。
+
+代码位置：`src/cryo_tank/solver.py:70-124`

src/cryo_tank/README.md

@@ -0,0 +1,381 @@
+# 低温液氮储箱仿真模块使用说明
+
+本目录实现了一个低温液氮储箱模型，用于模拟液氮储箱在氦气增压、液氮进出口流量、环境漏热和气液两区换热作用下的温度、液位、压力、氦气用量等随时间变化。
+
+模型入口文件是 `main.py`，主要参数集中在 `config.py`。运行后会把数据和图像写入 `results/cryo_tank/`。
+
+## 运行方法
+
+在项目根目录运行：
+
+```bash
+python3 src/cryo_tank/main.py
+```
+
+如果使用本地虚拟环境，先安装依赖：
+
+```bash
+.venv/bin/python -m pip install numpy scipy matplotlib pillow CoolProp pytest python-docx
+```
+
+然后运行：
+
+```bash
+.venv/bin/python src/cryo_tank/main.py
+```
+
+建议从项目根目录运行，而不是进入 `src/cryo_tank/` 后运行。这样输出路径会稳定写到：
+
+```text
+results/cryo_tank/
+```
+
+## 输出文件
+
+默认输出目录由 `config.py` 中的 `OUTPUT_DIR` 控制：
+
+```python
+OUTPUT_DIR = "results/cryo_tank"
+```
+
+运行完成后会生成：
+
+- `cryo_tank_history.npz`：压缩后的完整时序数据，适合用 Python 后处理。
+- `cryo_tank_history.csv`：CSV 表格数据，适合用 Excel、Origin 或其他工具查看。
+- `cryo_tank_temperatures.png`：液相温度和气枕区温度曲线。
+- `cryo_tank_level.png`：液位和充满率曲线。
+- `cryo_tank_he_flow.png`：氦气质量和氦气流量曲线。
+- `cryo_tank_heat.png`：漏热和气液换热曲线。
+- `cryo_tank_pressure.png`：总压和氦气气枕压力曲线。
+
+如果没有看到 CSV 文件，先确认本地代码里是否有：
+
+```bash
+grep -n "save_history_csv" src/cryo_tank/main.py src/cryo_tank/output.py
+```
+
+然后确认实际输出位置：
+
+```bash
+find results/cryo_tank -maxdepth 1 -type f -print
+```
+
+## 代码文件功能
+
+### `main.py`
+
+仿真入口。主要流程是：
+
+1. 从 `config.py` 读取几何、工况、初始条件和求解参数。
+2. 创建漏热模型 `MLIHeatLeak`。
+3. 创建储箱模型 `CryoTank`。
+4. 调用 `solver.run()` 进行 ODE 求解。
+5. 调用 `output.py` 保存 `.npz`、`.csv` 和图像。
+
+如果只是运行已有工况，一般只需要执行这个文件，不需要改动它。
+
+### `config.py`
+
+集中存放默认参数。常用调参基本都在这个文件里完成，包括：
+
+- 储箱几何尺寸
+- 工作压力
+- 初始温度和初始气枕率
+- 液氮入口、出口流量
+- 氦气入口温度
+- 气液换热系数
+- 环境温度
+- 仿真结束时间
+- 求解器误差容限
+- 输出目录
+
+### `tank_model.py`
+
+核心物理模型，定义 `CryoTank` 类。
+
+状态量为：
+
+```text
+y = [m_liq, U_liq, T_ull]
+```
+
+含义分别是：
+
+- `m_liq`：液氮质量
+- `U_liq`：液相内能
+- `T_ull`：气枕区温度
+
+`derive(y)` 会根据状态量计算温度、体积、液位、充满率、氦气质量、气枕区总内能和氦气气枕压力等派生量。当前模型按恒压纯氦气枕处理，`m_He = rho_He(T_ull, P_work) * V_ull`，所以氦气质量是派生量，不是独立 ODE 状态。
+
+`rhs(t, y)` 是 ODE 右端函数，由 `scipy.integrate.solve_ivp` 调用。
+
+### `solver.py`
+
+ODE 求解器封装。主要函数是：
+
+```python
+run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0)
+```
+
+它会：
+
+- 调用 `tank.initial_state()` 获取初始状态。
+- 使用 `solve_ivp` 求解。
+- 在液氮质量降到 0 时终止计算。
+- 对每个输出时刻计算温度、压力、液位、漏热、换热等时序数据。
+- 返回 `history` 字典。
+
+### `properties.py`
+
+物性计算模块。
+
+主要功能：
+
+- 液氮密度、焓、内能计算。
+- 根据液氮比内能和压力反推液相温度。
+- 氦气密度、焓、内能、比热计算。
+- 根据氦气压力和密度反推气枕温度。
+- 计算氦气恒压密度温度导数 `(partial rho / partial T)_P`。
+- 计算氦气恒压内能温度导数 `(partial u / partial T)_P`。
+
+该模块依赖 `CoolProp`。
+
+### `heat_leak.py`
+
+漏热模型模块。
+
+当前包含：
+
+- `HeatLeakModel`：漏热模型基类。
+- `MLIHeatLeak`：真空多层绝热模型，按固定热流密度 `q_mli` 计算。
+- `FoamHeatLeak`：泡沫或包覆绝热模型，按一维稳态导热计算。
+
+默认在 `main.py` 中使用：
+
+```python
+heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0)
+```
+
+### `output.py`
+
+输出模块。
+
+负责保存：
+
+- 压缩时序数据 `.npz`
+- 表格数据 `.csv`
+- 温度、液位、氦气流量、热流和压力图像
+
+绘图使用 `matplotlib` 的 `Agg` 后端，因此不需要图形界面。
+
+## 常用调参方法
+
+### 修改仿真时间
+
+在 `config.py` 中修改：
+
+```python
+T_END = 3600.0
+```
+
+单位是秒。比如仿真 10 分钟：
+
+```python
+T_END = 600.0
+```
+
+### 修改储箱尺寸
+
+在 `config.py` 中修改：
+
+```python
+V_TOTAL = 420.1e-3
+H_TANK = 0.5
+```
+
+其中：
+
+- `V_TOTAL` 是储箱总体积，单位 `m^3`。
+- `H_TANK` 是圆柱储箱高度，单位 `m`。
+
+横截面积、直径、侧面积和总面积会自动由这两个参数计算。
+
+### 修改工作压力
+
+在 `config.py` 中修改：
+
+```python
+P_WORKING = 0.17e6
+P_MAX = 0.8e6
+```
+
+单位是 Pa，代码中使用绝对压力。
+
+`P_WORKING` 是模型维持的工作压力，`P_MAX` 是最大承压限制，用于参数合法性检查。
+
+### 修改初始充液量
+
+在 `config.py` 中修改：
+
+```python
+ULLAGE_FRACTION = 0.30
+```
+
+该参数表示初始气枕体积分数。比如：
+
+- `0.30` 表示气枕区占 30%，液氮占 70%。
+- `0.20` 表示气枕区占 20%，液氮占 80%。
+
+取值必须在 0 到 1 之间。
+
+### 修改液氮进出口流量
+
+在 `config.py` 中修改：
+
+```python
+MDOT_IN_LN2 = 1.144
+MDOT_OUT_LN2 = 1.1895
+```
+
+单位是 `kg/s`。
+
+当前液氮质量变化率为：
+
+```text
+dm_liq/dt = MDOT_IN_LN2 - MDOT_OUT_LN2
+```
+
+如果出口流量大于入口流量，液位会逐渐下降。
+
+### 修改入口温度
+
+在 `config.py` 中修改：
+
+```python
+T_IN_LN2 = 77.0
+T_IN_HE = 100.0
+```
+
+单位是 K。
+
+### 修改环境漏热
+
+默认漏热模型在 `main.py` 中设置：
+
+```python
+heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0)
+```
+
+其中 `q_mli` 是单位面积漏热，单位 `W/m^2`。增大 `q_mli` 会增加外界进入储箱的热量。
+
+如果需要改用泡沫绝热模型，可以在 `main.py` 中引入并替换为：
+
+```python
+from cryo_tank.heat_leak import FoamHeatLeak
+
+heat_leak = FoamHeatLeak(A_total=A_TOTAL, k_eff=0.03, delta=0.05)
+```
+
+其中：
+
+- `k_eff` 是等效导热系数，单位 `W/(m*K)`。
+- `delta` 是绝热层厚度，单位 `m`。
+
+### 修改气液换热
+
+在 `config.py` 中修改：
+
+```python
+H_CONV_SURFACE = 50.0
+```
+
+单位是 `W/(m^2*K)`。该参数控制液相和气枕区之间的换热强度。
+
+### 修改求解精度
+
+在 `config.py` 中修改：
+
+```python
+RTOL = 1e-8
+ATOL = 1e-10
+```
+
+如果计算较慢，可以适当放宽；如果需要更高精度，可以适当收紧。收紧误差容限通常会增加计算时间。
+
+`solver.run()` 还有一个默认参数：
+
+```python
+max_step=10.0
+```
+
+表示最大时间步长，单位秒。如果需要更密集的输出点，可以在调用 `run()` 时减小它。
+
+## 后处理数据
+
+读取 `.npz` 数据示例：
+
+```python
+import numpy as np
+
+data = np.load("results/cryo_tank/cryo_tank_history.npz")
+print(data.files)
+print(data["t"])
+print(data["T_liq"])
+```
+
+读取 CSV 数据示例：
+
+```python
+import numpy as np
+
+data = np.loadtxt(
+    "results/cryo_tank/cryo_tank_history.csv",
+    delimiter=",",
+    skiprows=1,
+)
+print(data.shape)
+```
+
+CSV 第一行是列名，列名来自 `history` 字典，包括：
+
+- `t`
+- `m_liq`
+- `U_liq`
+- `U_ull`
+- `T_liq`
+- `T_ull`
+- `m_He`
+- `mdot_He`
+- `V_liq`
+- `V_ull`
+- `liquid_level`
+- `fill_fraction`
+- `P_He`
+- `P_total`
+- `Q_leak`
+- `Q_leak_liq`
+- `Q_leak_ull`
+- `Q_liq_to_ull`
+
+## 测试
+
+运行低温储箱相关测试：
+
+```bash
+pytest -q tests/cryo_tank
+```
+
+运行全项目测试：
+
+```bash
+pytest -q
+```
+
+## 注意事项
+
+- 本模型当前将储箱分为液相区和气枕区两个区域，不是完整 CFD 模型。
+- 气枕区按纯氦气处理，物性由 CoolProp 计算。
+- 气枕温度 `T_ull` 是 ODE 状态量；氦气质量 `m_He` 和补气流量 `mdot_He` 是恒压条件下的派生结果。
+- 液氮物性依赖 `CoolProp`，缺少该包会导致程序无法启动。
+- `results/` 目录默认被 `.gitignore` 忽略，仿真输出不会自动上传到 Git。
+- 建议每次修改参数后保存对应工况说明，避免不同仿真结果混在同一个输出目录中。

src/cryo_tank/config.py

@@ -3,14 +3,8 @@
 Configuration for the cryogenic LN2 tank simulation.
 Pure data module -- no functions, no side effects.
 """
-import math
 
-# ---------- Gas constants ----------
+import math
-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)
@@ -36,11 +30,11 @@ 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
+H_CONV_SURFACE = 0.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)
+T_END = 4818.884042  # s  (liquid level reaches 5%)
 RTOL = 1e-8
 ATOL = 1e-10
 

src/cryo_tank/main.py

@@ -21,7 +21,7 @@ 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,
+    save_history, save_history_csv, plot_temperatures, plot_liquid_level,
     plot_he_flow, plot_heat_fluxes, plot_pressure,
 )
 
@@ -49,7 +49,7 @@ def main():
     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"  P_He = {info0['P_He']/1e6:.4f} MPa")
     print(f"Running to t_end = {T_END:.0f} s ...")
     print()
 
@@ -67,6 +67,7 @@ def main():
 
     # --- Output ---
     save_history(history, os.path.join(OUTPUT_DIR, "cryo_tank_history.npz"))
+    save_history_csv(history, os.path.join(OUTPUT_DIR, "cryo_tank_history.csv"))
     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"))

src/cryo_tank/output.py

@@ -1,7 +1,7 @@
 # src/cryo_tank/output.py
 """
 Output helpers for the cryogenic tank simulation.
-Generates PNG plots and NPZ data files.
+Generates PNG plots, NPZ data files, and CSV tables.
 """
 import os
 
@@ -19,6 +19,27 @@ def save_history(history, path):
     np.savez_compressed(path, **history)
 
 
+def save_history_csv(history, path):
+    """Save all history time series to a CSV file."""
+    dirname = os.path.dirname(path)
+    if dirname:
+        os.makedirs(dirname, exist_ok=True)
+
+    field_names = list(history.keys())
+    columns = [np.asarray(history[name]) for name in field_names]
+    n_rows = len(columns[0]) if columns else 0
+
+    for name, values in zip(field_names, columns):
+        if values.ndim != 1:
+            raise ValueError(f"CSV history field must be 1D: {name}")
+        if len(values) != n_rows:
+            raise ValueError(f"CSV history field length mismatch: {name}")
+
+    data = np.column_stack(columns) if columns else np.empty((0, 0))
+    header = ",".join(field_names)
+    np.savetxt(path, data, delimiter=",", header=header, comments="")
+
+
 def plot_temperatures(history, path):
     """Plot T_liq and T_ull vs time."""
     fig, ax = plt.subplots(figsize=(10, 5))
@@ -95,12 +116,11 @@ def plot_heat_fluxes(history, path):
 
 
 def plot_pressure(history, path):
-    """Plot tank pressure (P_total, P_N2, P_He) vs time."""
+    """Plot tank pressure (P_total and 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]')

src/cryo_tank/properties.py

@@ -1,13 +1,11 @@
 # src/cryo_tank/properties.py
 """
-Fluid property wrappers for liquid nitrogen, N2 vapor, and helium.
+Fluid property wrappers for liquid nitrogen and helium.
 
-Performance strategy (per spec Section 8.1):
+Performance strategy:
-  - N2 (liquid & vapor): CoolProp with persistent AbstractState objects
+  - N2 liquid: CoolProp with a persistent AbstractState object
-  - He: analytical ideal gas (cp=5196.2 J/(kg*K), cv=3117.1 J/(kg*K))
+  - He: CoolProp with a persistent AbstractState object
-  - Lookup tables for the ODE hot path (built at import time)
 """
-import numpy as np
 import CoolProp.CoolProp as CP
 from CoolProp import AbstractState
 
@@ -41,81 +39,85 @@ def ln2_u(T, P):
 
 
 def ln2_T_from_u(u, P):
-    """Recover LN2 temperature from specific internal energy [K].
+    """Recover LN2 temperature from pressure and specific internal energy [K]."""
+    _n2_state.update(CP.PUmass_INPUTS, P, u)
+    return _n2_state.T()
 
-    Uses lookup table interpolation for speed; falls back to CoolProp
+
-    if outside the table range.
+def ln2_drho_dT_const_p(T, P):
-    """
+    """LN2 density temperature derivative at constant pressure [kg/(m^3*K)]."""
-    return float(np.interp(u, _ln2_u_table, _ln2_T_table))
+    _n2_state.update(CP.PT_INPUTS, P, T)
+    return _n2_state.first_partial_deriv(CP.iDmass, CP.iT, CP.iP)
+
+
+def ln2_du_dT_const_p(T, P):
+    """LN2 internal-energy temperature derivative at constant pressure [J/(kg*K)]."""
+    _n2_state.update(CP.PT_INPUTS, P, T)
+    return _n2_state.first_partial_deriv(CP.iUmass, CP.iT, CP.iP)
 
 
 # ---------------------------------------------------------------------------
-# N2 vapor properties (at saturation or specified conditions)
+# Helium properties via CoolProp
 # ---------------------------------------------------------------------------
-def n2_sat_pressure(T):
+_HE_P_REF = 170000.0
-    """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):
+def he_rho(T, P=_HE_P_REF):
-    """N2 saturated vapor specific internal energy [J/kg] at temperature T."""
+    """He density [kg/m^3]."""
-    _n2_state.update(CP.QT_INPUTS, 1.0, T)
+    _he_state.update(CP.PT_INPUTS, P, T)
-    return _n2_state.umass()
+    return _he_state.rhomass()
 
 
-def n2_vapor_rho(T):
+def he_cp(T=100.0, P=_HE_P_REF):
-    """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
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.cpmass()
 
 
-def he_cv():
+def he_cv(T=100.0, P=_HE_P_REF):
     """He specific heat at constant volume [J/(kg*K)]."""
-    return _HE_CV
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.cvmass()
 
 
-def he_h(T):
+def he_h(T, P=_HE_P_REF):
-    """He specific enthalpy [J/kg] (ideal gas)."""
+    """He specific enthalpy [J/kg]."""
-    return _HE_CP * T + _HE_H_REF
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.hmass()
 
 
-def he_u(T):
+def he_u(T, P=_HE_P_REF):
-    """He specific internal energy [J/kg] (ideal gas)."""
+    """He specific internal energy [J/kg]."""
-    return _HE_CV * T + _HE_U_REF
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.umass()
 
 
-def he_T_from_u(u):
+def he_T_from_u(u, P=_HE_P_REF):
-    """Recover He temperature from specific internal energy [K]."""
+    """Recover He temperature from pressure and specific internal energy [K]."""
-    return (u - _HE_U_REF) / _HE_CV
+    _he_state.update(CP.PUmass_INPUTS, P, u)
+    return _he_state.T()
 
 
-# ---------------------------------------------------------------------------
+def he_T_from_rho(P, rho):
-# Lookup table for LN2: u(T) -> T at P = 0.17 MPa (built at import time)
+    """Recover He temperature from pressure and density [K]."""
-# ---------------------------------------------------------------------------
+    _he_state.update(CP.DmassP_INPUTS, rho, P)
-_LN2_T_MIN = 65.0
+    return _he_state.T()
-_LN2_T_MAX = 82.0   # stay below saturation at 0.17 MPa (~82.03 K)
+
-_LN2_TABLE_N = 200
+
-_P_WORK = 170000.0
+def he_u_from_rho(P, rho):
+    """Recover He specific internal energy from pressure and density [J/kg]."""
+    _he_state.update(CP.DmassP_INPUTS, rho, P)
+    return _he_state.umass()
+
+
+def he_drho_dT_const_p(T, P=_HE_P_REF):
+    """He density temperature derivative at constant pressure [kg/(m^3*K)]."""
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.first_partial_deriv(CP.iDmass, CP.iT, CP.iP)
+
+
+def he_du_dT_const_p(T, P=_HE_P_REF):
+    """He internal-energy temperature derivative at constant pressure [J/(kg*K)]."""
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.first_partial_deriv(CP.iUmass, CP.iT, CP.iP)
 
-_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

@@ -71,7 +71,7 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
         't': t,
         'm_liq': sol.y[0],
         'U_liq': sol.y[1],
-        'U_ull': sol.y[2],
+        'U_ull': np.zeros(n),
         'T_liq': np.zeros(n),
         'T_ull': np.zeros(n),
         'm_He': np.zeros(n),
@@ -80,7 +80,6 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
         '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),
@@ -95,14 +94,14 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
 
         history['T_liq'][i] = info['T_liq']
         history['T_ull'][i] = info['T_ull']
+        history['U_ull'][i] = info['U_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']
+        history['P_total'][i] = info['P_He']
 
         # Recompute heat terms for recording
         T_liq = info['T_liq']
@@ -120,11 +119,9 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
         history['Q_leak'][i] = Q_leak
         history['Q_leak_liq'][i] = Q_leak_liq
         history['Q_leak_ull'][i] = Q_leak_ull
+        history['mdot_He'][i] = tank._solve_he_flow_rate(
+            info, Q_liq_to_ull, Q_leak_ull, Q_leak_liq
+        )
 
-    # 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

@@ -2,7 +2,7 @@
 """
 CryoTank: two-zone (liquid + ullage) cryogenic tank model.
 
-State vector y = [m_liq, U_liq, U_ull] (3 components).
+State vector y = [m_liq, U_liq, T_ull] (3 components).
 Derived quantities (T, V, m_He, etc.) computed by derive(y).
 ODE right-hand side provided by rhs(t, y).
 """
@@ -11,7 +11,6 @@ import warnings
 
 import numpy as np
 
-from cryo_tank.config import R_HE, R_N2
 from cryo_tank import properties as prop
 
 
@@ -60,7 +59,7 @@ class CryoTank:
 
         # 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)
+        self.h_in_he = prop.he_h(T_in_he, P_work)
 
         # Net liquid flow (constant)
         self.dm_liq_dt = mdot_in_ln2 - mdot_out_ln2
@@ -77,25 +76,17 @@ class CryoTank:
         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
+        # Ullage initial state: helium pressurization only.
-        # Use ideal gas law for N2 mass (consistent with derive/rhs which
+        # Nitrogen evaporation is intentionally not modeled, so the ullage
-        # treat ullage N2 as ideal gas: P_N2 = m_N2 * R_N2 * T / V)
+        # pressure is provided entirely by helium and the liquid sees the same
-        P_N2_0 = prop.n2_sat_pressure(T_init)
+        # tank pressure for property lookup.
-        self.m_N2_ull = P_N2_0 * V_ull_0 / (R_N2 * T_init)  # FIXED for all time
+        self._T_ull_0 = T_init
-
+        self._m_He_0 = prop.he_rho(T_init, P_work) * V_ull_0
-        P_He_0 = P_work - P_N2_0
+        self._U_ull_0 = self._m_He_0 * prop.he_u(T_init, P_work)
-        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 the ODE initial state vector y0 = [m_liq, U_liq, T_ull]."""
-        return np.array([self._m_liq_0, self._U_liq_0, self._U_ull_0])
+        return np.array([self._m_liq_0, self._U_liq_0, self._T_ull_0])
 
     def wetted_areas(self, liquid_level):
         """Return (A_wet, A_dry) for the given liquid level [m]."""
@@ -108,94 +99,45 @@ class CryoTank:
         """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.
+        fill_fraction, m_He, U_ull, P_He, etc. Helium provides the full
+        ullage pressure in this no-evaporation model.
         """
-        m_liq, U_liq, U_ull = y[0], y[1], y[2]
+        m_liq, U_liq, T_ull = y[0], y[1], y[2]
 
         # Liquid zone
-        u_liq = U_liq / m_liq  # specific internal energy
+        u_liq = U_liq / m_liq
         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
+        # Ullage zone: pure helium at tank pressure. T_ull is the ODE state;
+        # m_He is the amount required to maintain P_work in the current volume.
         V_ull = self.V_total - V_liq
-
+        P_He = self.P_work
-        # Solve T_ull from ullage energy
+        rho_He = prop.he_rho(T_ull, P_He)
-        T_ull = self._solve_ullage_temperature(U_ull, V_ull)
+        m_He = rho_He * V_ull
-
+        U_ull = m_He * prop.he_u(T_ull, P_He)
-        # 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 {
+            'm_liq': m_liq, 'U_liq': U_liq, 'u_liq': u_liq,
             '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,
+            'rho_liq': rho_liq, 'rho_He': rho_He,
-            'm_He': m_He, 'P_N2': P_N2, 'P_He': P_He,
+            'm_He': m_He, 'U_ull': U_ull, '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].
+        """ODE right-hand side: dy/dt = [dm_liq/dt, dU_liq/dt, dT_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)
@@ -206,77 +148,110 @@ class CryoTank:
         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._liquid_energy_rate(
-        dU_liq_dt = (self.mdot_in_ln2 * self.h_in_ln2
+            info, Q_liq_to_ull, Q_leak_liq
+        )
+
+        # --- Ullage zone ---
+        dT_ull_dt, _ = self._solve_ullage_temperature_rate(
+            info, Q_liq_to_ull, Q_leak_ull, dU_liq_dt
+        )
+
+        return np.array([dm_liq_dt, dU_liq_dt, dT_ull_dt])
+
+    def _liquid_energy_rate(self, info, Q_liq_to_ull, Q_leak_liq):
+        """Return liquid-zone total internal-energy rate [W]."""
+        h_liq = prop.ln2_h(info['T_liq'], self.P_work)
+        return (self.mdot_in_ln2 * self.h_in_ln2
                 - self.mdot_out_ln2 * h_liq
                 - Q_liq_to_ull
                 + Q_leak_liq)
 
-        # --- Ullage zone ---
+    def _liquid_temperature_rate(self, info, dm_liq_dt, dU_liq_dt):
-        dU_ull_dt = mdot_He * self.h_in_he + Q_liq_to_ull + Q_leak_ull
+        """Return dT_liq/dt from u_liq=U_liq/m_liq at constant tank pressure."""
+        m_liq = info['m_liq']
+        u_liq = info['u_liq']
+        du_dT = prop.ln2_du_dT_const_p(info['T_liq'], self.P_work)
+        denominator = m_liq * du_dT
+        if abs(denominator) < 1e-30:
+            return 0.0
+        return (dU_liq_dt - u_liq * dm_liq_dt) / denominator
 
-        # --- Warnings ---
+    def _ullage_volume_rate(self, info, dm_liq_dt, dU_liq_dt):
-        if T_liq > self._T_sat - 1.0:
+        """Return dV_ull/dt including liquid density variation with temperature."""
-            warnings.warn(
+        rho_liq = info['rho_liq']
-                f"T_liq={T_liq:.2f}K approaching saturation ({self._T_sat:.1f}K); "
+        m_liq = info['m_liq']
-                "evaporation effects may be significant."
+        T_liq = info['T_liq']
+        dT_liq_dt = self._liquid_temperature_rate(
+            info, dm_liq_dt, dU_liq_dt
         )
+        drho_dT = prop.ln2_drho_dT_const_p(T_liq, self.P_work)
 
-        return np.array([dm_liq_dt, dU_liq_dt, dU_ull_dt])
+        dV_liq_dt = (dm_liq_dt / rho_liq
+                     - m_liq * drho_dT * dT_liq_dt / rho_liq ** 2)
+        return -dV_liq_dt
 
-    def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull):
+    def _he_mass_energy_partials(self, info):
-        """Analytically solve for m_dot_He from dP/dt = 0 constraint.
+        """Return local partials for m(T,V) and U(T,V) at constant pressure."""
-
-        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 = info['rho_He']
-        rho_liq = info['rho_liq']
+        u = prop.he_u(T_ull, self.P_work)
+        drho_dT = prop.he_drho_dT_const_p(T_ull, self.P_work)
+        du_dT = prop.he_du_dT_const_p(T_ull, self.P_work)
 
-        # dV_ull/dt = -dV_liq/dt = -dm_liq/dt / rho_liq = -(mdot_in - mdot_out) / rho_liq
+        m_T = V_ull * drho_dT
-        dV_ull_dt = -self.dm_liq_dt / rho_liq  # positive when liquid drains
+        m_V = rho
+        U_T = V_ull * (u * drho_dT + rho * du_dT)
+        U_V = rho * u
+        return m_T, U_T, m_V, U_V
 
-        # Total ullage heat input (excluding He inlet, which we're solving for)
+    def _solve_ullage_temperature_rate(self, info, Q_liq_to_ull, Q_leak_ull,
+                                       dU_liq_dt):
+        """Solve dT_ull/dt and He inlet flow from constant-pressure EOS.
+
+        T_ull is the ODE state. The pure-He ullage is maintained at P_work,
+        so m_He = rho(P_work, T_ull) * V_ull is a derived quantity.
+        """
+        dV_ull_dt = self._ullage_volume_rate(
+            info, self.dm_liq_dt, dU_liq_dt
+        )
         Q_ull_no_he = Q_liq_to_ull + Q_leak_ull
 
-        # Ullage total cv*mass (for dT_ull/dt estimation)
+        m_T, U_T, m_V, U_V = self._he_mass_energy_partials(info)
-        cv_N2 = 743.0  # J/(kg*K), N2 vapor
+        denominator = U_T - self.h_in_he * m_T
-        cv_He = prop.he_cv()
+        if abs(denominator) < 1e-30:
-        C_ull = self.m_N2_ull * cv_N2 + m_He * cv_He  # total heat capacity [J/K]
+            return 0.0, 0.0
 
-        R_mix = self.m_N2_ull * R_N2 + m_He * R_HE  # effective "mR" [J/K]
+        dT_ull_dt = (Q_ull_no_he
+                     - (U_V + self.P_work - self.h_in_he * m_V)
+                     * dV_ull_dt) / denominator
+        mdot_He = m_T * dT_ull_dt + m_V * dV_ull_dt
 
-        # 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."
+                "clamping reported flow to 0. Pressure control may be invalid."
             )
             mdot_He = 0.0
 
+        return dT_ull_dt, mdot_He
+
+    def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull,
+                            Q_leak_liq=None):
+        """Return derived He inlet mass flow for the current state [kg/s]."""
+        if Q_leak_liq is None:
+            Q_leak = self.heat_leak_model.compute(info['T_liq'], self.T_env)
+            A_wet, A_dry = self.wetted_areas(info['liquid_level'])
+            A_total = A_wet + A_dry
+            Q_leak_liq = Q_leak * A_wet / A_total if A_total > 0 else 0.0
+
+        dU_liq_dt = self._liquid_energy_rate(
+            info, Q_liq_to_ull, Q_leak_liq
+        )
+        _, mdot_He = self._solve_ullage_temperature_rate(
+            info, Q_liq_to_ull, Q_leak_ull, dU_liq_dt
+        )
         return mdot_He

src/cylinder/Index.md

@@ -0,0 +1,338 @@
+## 1. 模型定位
+
+`gas_cylinder.py` 实现的是高压气瓶的 0D 集总参数模型 `HighPressureGasCylinder`。模型只描述气瓶内部平均热力状态，不解析气瓶内的轴向、径向或局部流场分布。
+
+主要假设：
+
+- 气瓶有效容积 `V` 固定。
+- 气瓶内工质均匀，使用单一平均密度和单一平均比内能描述。
+- 独立守恒状态只有总质量 `mass` 和总内能 `U`。
+- 压力、温度和比焓通过 CoolProp 由当前守恒状态恢复。
+- 模型不自行计算边界流量；外部模型给出质量流率和能量流率后，气瓶只更新守恒量。
+
+对应代码：`src/cylinder/gas_cylinder.py:20-142`
+
+## 2. 参数与状态
+
+默认参数定义在 `src/cylinder/gas_cylinder.py:14-17`：
+
+| 参数 | 默认值 | 含义 |
+| --- | ---: | --- |
+| `DEFAULT_FLUID` | `Helium` | CoolProp 工质名称 |
+| `DEFAULT_P_INIT` | `18.031e6 Pa` | 初始绝对压力 |
+| `DEFAULT_VOLUME` | `20e-3 m^3` | 气瓶容积，20 L |
+| `DEFAULT_T_INIT` | `78.0 K` | 初始温度 |
+
+构造函数输入为：
+
+$$
+\left(\mathrm{fluid},\ P_0,\ V,\ T_0,\ \mathrm{backend}\right)
+$$
+
+其中 `backend` 默认是 `HEOS`。构造函数会检查：
+
+$$
+\mathrm{fluid}\ne\varnothing,\quad P_0>0,\quad V>0,\quad T_0>0
+$$
+
+气瓶模型的独立状态向量为：
+
+$$
+\mathbf{y}_{\mathrm{cyl}}=\left[m_{\mathrm{cyl}},\ U_{\mathrm{cyl}}\right]
+$$
+
+其中：
+
+- `mass` 对应 $m_{\mathrm{cyl}}$，单位 kg。
+- `U` 对应 $U_{\mathrm{cyl}}$，单位 J。
+- 压力 `P`、温度 `T`、密度 `rho`、比内能 `specific_internal_energy` 和比焓 `h` 都是由这两个守恒量派生出来的量。
+
+## 3. 初始状态计算
+
+初始化发生在 `src/cylinder/gas_cylinder.py:51-56`。计算流程如下。
+
+第一步，建立 CoolProp 状态对象：
+
+```python
+self._state = AbstractState(backend, fluid)
+```
+
+第二步，用初始压力和初始温度更新热力状态：
+
+```python
+self._state.update(CP.PT_INPUTS, P_init, T_init)
+```
+
+对应数学表达为：
+
+$$
+\mathrm{state}_0=\mathrm{CoolProp}\left(P_0,T_0,\mathrm{fluid},\mathrm{backend}\right)
+$$
+
+第三步，从 CoolProp 读取初始密度：
+
+$$
+\rho_0=\rho\left(P_0,T_0\right)
+$$
+
+代码中对应：
+
+```python
+rho = self._state.rhomass()
+```
+
+第四步，由固定容积计算气瓶初始总质量：
+
+$$
+m_{\mathrm{cyl},0}=\rho_0 V
+$$
+
+代码中对应：
+
+```python
+self.mass = rho * V
+```
+
+第五步，从 CoolProp 读取初始比内能，并计算总内能：
+
+$$
+U_{\mathrm{cyl},0}=m_{\mathrm{cyl},0}u_0
+$$
+
+其中：
+
+$$
+u_0=u\left(P_0,T_0\right)
+$$
+
+代码中对应：
+
+```python
+self.U = self.mass * self._state.umass()
+```
+
+因此，模型初始化后不再单独保存初始压力和初始温度；后续压力和温度都通过当前 `mass`、`U` 和 `V` 重新计算。
+
+## 4. 状态恢复计算
+
+运行过程中，模型每次需要热力性质时都会调用 `_update_state()`，代码位置为 `src/cylinder/gas_cylinder.py:58-62`。
+
+### 4.1 密度
+
+气瓶体积固定，因此当前平均密度为：
+
+$$
+\rho_{\mathrm{cyl}}=\frac{m_{\mathrm{cyl}}}{V}
+$$
+
+代码位置：`src/cylinder/gas_cylinder.py:64-67`
+
+### 4.2 比内能
+
+当前比内能由总内能除以总质量得到：
+
+$$
+u_{\mathrm{cyl}}=\frac{U_{\mathrm{cyl}}}{m_{\mathrm{cyl}}}
+$$
+
+代码位置：`src/cylinder/gas_cylinder.py:69-72`
+
+这一步要求 `mass > 0`。`apply_flux()` 在每次通量更新后会检查质量不能为非正值。
+
+### 4.3 CoolProp 输入对
+
+模型用质量密度和质量比内能作为 CoolProp 输入对：
+
+```python
+self._state.update(CP.DmassUmass_INPUTS, rho, u)
+```
+
+对应数学表达为：
+
+$$
+\mathrm{state}=\mathrm{CoolProp}\left(\rho_{\mathrm{cyl}},u_{\mathrm{cyl}}\right)
+$$
+
+这里的关键点是：`P` 和 `T` 不是独立积分变量，而是 CoolProp 根据 $(\rho,u)$ 反算出来的热力结果。
+
+### 4.4 压力、温度和比焓
+
+状态恢复后，模型提供三个主要派生量：
+
+$$
+\begin{aligned}
+P_{\mathrm{cyl}} &= P\left(\rho_{\mathrm{cyl}},u_{\mathrm{cyl}}\right) \\
+T_{\mathrm{cyl}} &= T\left(\rho_{\mathrm{cyl}},u_{\mathrm{cyl}}\right) \\
+h_{\mathrm{cyl}} &= h\left(\rho_{\mathrm{cyl}},u_{\mathrm{cyl}}\right)
+\end{aligned}
+$$
+
+对应代码：
+
+- `P`：`src/cylinder/gas_cylinder.py:74-77`
+- `T`：`src/cylinder/gas_cylinder.py:79-82`
+- `h`：`src/cylinder/gas_cylinder.py:84-87`
+
+当前模型只保留上述三个热力派生量，不再保留额外的热物性辅助输出。
+
+## 5. 边界 ghost state 计算
+
+`ghost_state()` 用于给一维守恒流动模型提供气瓶侧边界状态，代码位置为 `src/cylinder/gas_cylinder.py:89-97`。
+
+返回向量为：
+
+$$
+\mathbf{q}_{\mathrm{ghost}}=
+\left[
+\rho_{\mathrm{cyl}},\
+0,\
+\rho_{\mathrm{cyl}}u_{\mathrm{cyl}}
+\right]^T
+$$
+
+三个分量分别是：
+
+| 分量 | 代码表达 | 物理意义 |
+| --- | --- | --- |
+| $q_1$ | `self.rho` | 质量密度 |
+| $q_2$ | `0.0` | 动量密度 |
+| $q_3$ | `self.rho * self.specific_internal_energy` | 总能量密度 |
+
+由于气瓶是 0D 静止控制体，边界 ghost state 的速度取 0，所以动量密度为 0。总能量密度中也没有宏观动能项：
+
+$$
+\rho E = \rho u + \frac{1}{2}\rho v^2 = \rho u \quad \left(v=0\right)
+$$
+
+注意：这个 ghost state 使用的是真实流体内能。如果下游管路求解器使用理想气体状态方程，耦合前需要确认状态方程和能量定义是否一致。
+
+## 6. 边界通量更新计算
+
+`apply_flux(mdot, edot, dt, sign)` 根据外部边界通量更新气瓶守恒量，代码位置为 `src/cylinder/gas_cylinder.py:99-128`。
+
+输入量含义：
+
+| 参数 | 单位 | 含义 |
+| --- | --- | --- |
+| `mdot` | kg/s | 边界质量流率 |
+| `edot` | W | 边界能量流率 |
+| `dt` | s | 时间步长 |
+| `sign` | - | 通量方向符号 |
+
+方向符号定义为：
+
+$$
+s=
+\begin{cases}
+-1, & \text{气体从气瓶流出} \\
++1, & \text{气体流入气瓶}
+\end{cases}
+$$
+
+离散更新方程为：
+
+$$
+\begin{aligned}
+m_{\mathrm{cyl}}^{n+1} &= m_{\mathrm{cyl}}^{n}+s\dot{m}\Delta t \\
+U_{\mathrm{cyl}}^{n+1} &= U_{\mathrm{cyl}}^{n}+s\dot{E}\Delta t
+\end{aligned}
+$$
+
+代码中对应：
+
+```python
+self.mass += sign * mdot * dt
+self.U += sign * edot * dt
+```
+
+连续形式为：
+
+$$
+\begin{aligned}
+\frac{dm_{\mathrm{cyl}}}{dt} &= s\dot{m} \\
+\frac{dU_{\mathrm{cyl}}}{dt} &= s\dot{E}
+\end{aligned}
+$$
+
+如果外部模型使用边界比焓计算能量流率，常见形式为：
+
+$$
+\dot{E}=\dot{m}h_{\mathrm{boundary}}
+$$
+
+但气瓶模型本身不强制 `edot = mdot * h`，只接收外部传入的能量流率。这使它可以用于不同边界模型，但也要求调用方保证质量通量和能量通量的物理一致性。
+
+### 6.1 数值保护
+
+通量更新前后包含三类保护：
+
+$$
+\Delta t\ge0,\quad s\in\{-1,+1\},\quad m_{\mathrm{cyl}}^{n+1}>0
+$$
+
+对应行为：
+
+- `dt < 0` 时抛出 `ValueError`。
+- `sign` 不是 `-1` 或 `+1` 时抛出 `ValueError`。
+- 更新后 `mass <= 0` 时抛出 `RuntimeError`。
+- 更新完成后立即调用 `_update_state()`，让 CoolProp 立刻检查新的 $(\rho,u)$ 是否处于有效热力状态。
+
+## 7. 状态汇总
+
+`state_summary()` 返回当前气瓶状态字典，代码位置为 `src/cylinder/gas_cylinder.py:130-142`。
+
+返回字段为：
+
+| 字段 | 含义 |
+| --- | --- |
+| `fluid` | 工质名称 |
+| `V` | 气瓶固定容积 |
+| `mass` | 当前总质量 |
+| `rho` | 当前平均密度 |
+| `P` | 当前绝对压力 |
+| `T` | 当前温度 |
+| `U` | 当前总内能 |
+| `u` | 当前比内能 |
+| `h` | 当前比焓 |
+
+该函数主要用于后处理、日志输出或快速检查模型状态。
+
+## 8. 耦合使用简述
+
+在低温贮箱和气瓶耦合示例中，气瓶状态通常作为整体 ODE 状态向量的一部分：
+
+$$
+\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}},\ m_{\mathrm{cyl}},\ U_{\mathrm{cyl}}\right]
+$$
+
+气瓶向贮箱供气时，可写为：
+
+$$
+\frac{dm_{\mathrm{cyl}}}{dt}=-\dot{m}_{\mathrm{to\ tank}},\quad
+\frac{dU_{\mathrm{cyl}}}{dt}=-\dot{E}_{\mathrm{to\ tank}}
+$$
+
+并可用气瓶压力相对贮箱工作压力的差值作为终止条件：
+
+$$
+g(t,\mathbf{y})=P_{\mathrm{cyl}}-P_{\mathrm{work}}
+$$
+
+示例代码位置：`examples/cryo_tank_cylinder_system.py`
+
+## 9. 运行和测试
+
+气瓶模型可直接导入使用：
+
+```python
+from cylinder import HighPressureGasCylinder
+
+cylinder = HighPressureGasCylinder()
+summary = cylinder.state_summary()
+```
+
+相关测试覆盖默认状态、可配置输入、状态汇总字段、ghost state、出流通量更新和非法输入校验。测试命令：
+
+```bash
+pytest -q tests/test_gas_cylinder.py
+```

src/cylinder/__init__.py

@@ -0,0 +1,17 @@
+"""High-pressure gas cylinder models."""
+
+from cylinder.gas_cylinder import (
+    DEFAULT_FLUID,
+    DEFAULT_P_INIT,
+    DEFAULT_T_INIT,
+    DEFAULT_VOLUME,
+    HighPressureGasCylinder,
+)
+
+__all__ = [
+    "DEFAULT_FLUID",
+    "DEFAULT_P_INIT",
+    "DEFAULT_T_INIT",
+    "DEFAULT_VOLUME",
+    "HighPressureGasCylinder",
+]

src/cylinder/gas_cylinder.py

@@ -0,0 +1,142 @@
+"""
+High-pressure gas cylinder model with configurable working fluid.
+
+The cylinder is a 0D lumped model. Its primary state is total mass and
+total internal energy; pressure, temperature, density, and enthalpy are
+recovered from CoolProp on demand.
+"""
+
+import numpy as np
+import CoolProp.CoolProp as CP
+from CoolProp import AbstractState
+
+
+DEFAULT_FLUID = "Helium"
+DEFAULT_P_INIT = 18.031e6   # Pa, absolute
+DEFAULT_VOLUME = 20e-3      # m^3, 20 L
+DEFAULT_T_INIT = 78.0       # K
+
+
+class HighPressureGasCylinder:
+    """0D high-pressure gas cylinder.
+
+    Parameters
+    ----------
+    fluid : str
+        CoolProp fluid name, for example ``"Helium"`` or ``"Nitrogen"``.
+    P_init : float
+        Initial absolute pressure [Pa].
+    V : float
+        Cylinder volume [m^3].
+    T_init : float
+        Initial temperature [K].
+    backend : str
+        CoolProp backend name. ``"HEOS"`` is used by default.
+    """
+
+    def __init__(self, fluid=DEFAULT_FLUID, P_init=DEFAULT_P_INIT,
+                 V=DEFAULT_VOLUME, T_init=DEFAULT_T_INIT, backend="HEOS"):
+        if not fluid:
+            raise ValueError("fluid must be a non-empty CoolProp fluid name")
+        if P_init <= 0:
+            raise ValueError("P_init must be > 0")
+        if V <= 0:
+            raise ValueError("V must be > 0")
+        if T_init <= 0:
+            raise ValueError("T_init must be > 0")
+
+        self.fluid = fluid
+        self.backend = backend
+        self.V = V
+        self._state = AbstractState(backend, fluid)
+
+        self._state.update(CP.PT_INPUTS, P_init, T_init)
+        rho = self._state.rhomass()
+        self.mass = rho * V
+        self.U = self.mass * self._state.umass()
+
+    def _update_state(self):
+        rho = self.rho
+        u = self.specific_internal_energy
+        self._state.update(CP.DmassUmass_INPUTS, rho, u)
+        return self._state
+
+    @property
+    def rho(self):
+        """Gas density [kg/m^3]."""
+        return self.mass / self.V
+
+    @property
+    def specific_internal_energy(self):
+        """Specific internal energy [J/kg]."""
+        return self.U / self.mass
+
+    @property
+    def P(self):
+        """Absolute pressure [Pa]."""
+        return self._update_state().p()
+
+    @property
+    def T(self):
+        """Temperature [K]."""
+        return self._update_state().T()
+
+    @property
+    def h(self):
+        """Specific enthalpy [J/kg]."""
+        return self._update_state().hmass()
+
+    def ghost_state(self):
+        """Return a stagnant conservative state ``[rho, rho*u, rho*E]``.
+
+        This is useful as a boundary ghost cell for conservative flow models.
+        The state uses the real-fluid internal energy density. A pipe solver
+        that assumes an ideal-gas EOS must still be checked for EOS consistency
+        before coupling it directly to this real-fluid cylinder.
+        """
+        return np.array([self.rho, 0.0, self.rho * self.specific_internal_energy])
+
+    def apply_flux(self, mdot, edot, dt, sign):
+        """Update cylinder mass and energy from a boundary flux.
+
+        Parameters
+        ----------
+        mdot : float
+            Mass flow rate [kg/s].
+        edot : float
+            Energy flow rate [W].
+        dt : float
+            Time step [s].
+        sign : int
+            ``-1`` for outflow from the cylinder, ``+1`` for inflow.
+        """
+        if dt < 0:
+            raise ValueError("dt must be >= 0")
+        if sign not in (-1, 1):
+            raise ValueError("sign must be +1 or -1")
+
+        self.mass += sign * mdot * dt
+        self.U += sign * edot * dt
+
+        if self.mass <= 0:
+            raise RuntimeError(
+                f"Cylinder mass non-positive after apply_flux: mass={self.mass}, "
+                f"mdot={mdot}, edot={edot}, dt={dt}, sign={sign}"
+            )
+
+        # Force a thermodynamic validity check immediately after the update.
+        self._update_state()
+
+    def state_summary(self):
+        """Return common cylinder state quantities as a dictionary."""
+        return {
+            "fluid": self.fluid,
+            "V": self.V,
+            "mass": self.mass,
+            "rho": self.rho,
+            "P": self.P,
+            "T": self.T,
+            "U": self.U,
+            "u": self.specific_internal_energy,
+            "h": self.h,
+        }

src/main.py

@@ -1,130 +1,6 @@
-# src/main.py
+"""Compatibility entry point for the tank-pipe simulation."""
-"""
-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:
+from tank_pipe.main import main
-    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__":

src/tank_pipe/__init__.py

@@ -0,0 +1 @@
+"""0D-1D tank-pipe blowdown simulation package."""

src/tank_pipe/main.py

@@ -0,0 +1,130 @@
+# src/tank_pipe/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/tank_pipe/main.py
+"""
+import os
+import sys
+
+_HERE = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.dirname(_HERE))
+
+import numpy as np
+
+from tank_pipe.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_pipe.tank import Tank
+from tank_pipe.pipe import Pipe
+from tank_pipe.solver import run
+from tank_pipe.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/tank_pipe/output.py

@@ -1,4 +1,4 @@
-# src/output.py
+# src/tank_pipe/output.py
 """
 Output helpers: persistence (.npz), static plots (.png/.html), animation (.gif).
 Uses matplotlib's Agg backend so it works in headless environments.

src/tank_pipe/pipe.py

@@ -1,4 +1,4 @@
-# src/pipe.py
+# src/tank_pipe/pipe.py
 """
 1D finite-volume pipe for compressible Euler equations:
   dW/dt + dF(W)/dx = 0
@@ -12,8 +12,8 @@ Discretization:
 """
 
 import numpy as np
-from riemann import hll_flux, get_riemann_solver
+from tank_pipe.riemann import hll_flux, get_riemann_solver
-from friction import darcy_friction_factor
+from tank_pipe.friction import darcy_friction_factor
 
 
 class Pipe:

src/tank_pipe/riemann.py

@@ -1,4 +1,4 @@
-# src/riemann.py
+# src/tank_pipe/riemann.py
 """
 Riemann flux solvers for the 1D compressible Euler equations.
 

src/tank_pipe/tank.py

@@ -1,4 +1,4 @@
-# src/tank.py
+# src/tank_pipe/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,

tests/cryo_tank/test_properties.py

@@ -6,7 +6,7 @@ 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,
+    ln2_drho_dT_const_p, ln2_du_dT_const_p,
     he_u, he_h, he_cp, he_cv,
 )
 from cryo_tank.config import P_WORKING
@@ -33,25 +33,28 @@ class TestLN2Properties:
         T_recovered = ln2_T_from_u(u, P_WORKING)
         assert abs(T_recovered - T_orig) < 0.01
 
+    def test_ln2_density_derivative_matches_finite_difference(self):
+        T = 78.0
+        step = 1e-3
+        actual = ln2_drho_dT_const_p(T, P_WORKING)
+        expected = (
+            ln2_rho(T + step, P_WORKING)
+            - ln2_rho(T - step, P_WORKING)
+        ) / (2.0 * step)
+        assert actual < 0.0
+        assert abs(actual - expected) / abs(expected) < 1e-5
 
-class TestN2VaporProperties:
+    def test_ln2_internal_energy_derivative_is_positive(self):
-    """N2 vapor properties."""
+        du_dT = ln2_du_dT_const_p(78.0, P_WORKING)
-
+        assert du_dT > 0.0
-    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."""
+    """Helium properties."""
 
     def test_he_cp_near_5196(self):
         cp = he_cp()
-        assert abs(cp - 5196.2) < 10  # monatomic ideal gas
+        assert abs(cp - 5196.2) < 10
 
     def test_he_cv_near_3117(self):
         cv = he_cv()

tests/cryo_tank/test_tank_model.py

@@ -12,10 +12,9 @@ from cryo_tank.config import (
 )
 
 
-def _make_tank():
+def _make_tank(**overrides):
     """Create a CryoTank with default config and MLI heat leak."""
-    heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0)
+    kw = dict(
-    return CryoTank(
         V_total=V_TOTAL, H_tank=H_TANK,
         P_work=P_WORKING,
         T_init=T_INIT, ullage_fraction=ULLAGE_FRACTION,
@@ -23,8 +22,10 @@ def _make_tank():
         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,
+        heat_leak_model=MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0),
     )
+    kw.update(overrides)
+    return CryoTank(**kw)
 
 
 class TestGeometry:
@@ -69,10 +70,43 @@ class TestInitialState:
         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):
+    def test_initial_pressure_is_provided_entirely_by_helium(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
+        assert abs(P_He - P_WORKING) / P_WORKING < 1e-12
+
+
+class TestVolumeRate:
+
+    def test_ullage_volume_rate_includes_liquid_thermal_expansion(self):
+        tank = _make_tank(
+            mdot_in_ln2=0.0,
+            mdot_out_ln2=0.0,
+            h_conv=0.0,
+        )
+        y0 = tank.initial_state()
+        info = tank.derive(y0)
+
+        Q_liq_to_ull = 0.0
+        Q_leak = tank.heat_leak_model.compute(info['T_liq'], tank.T_env)
+        A_wet, A_dry = tank.wetted_areas(info['liquid_level'])
+        Q_leak_liq = Q_leak * A_wet / (A_wet + A_dry)
+        dU_liq_dt = tank._liquid_energy_rate(
+            info, Q_liq_to_ull, Q_leak_liq
+        )
+
+        dV_ull_dt = tank._ullage_volume_rate(
+            info, tank.dm_liq_dt, dU_liq_dt
+        )
+
+        y_next = y0.copy()
+        dt = 1.0
+        y_next[1] += dU_liq_dt * dt
+        finite_difference = (
+            tank.derive(y_next)['V_ull'] - info['V_ull']
+        ) / dt
+
+        assert dV_ull_dt < 0.0
+        assert abs(dV_ull_dt - finite_difference) / abs(finite_difference) < 1e-5

tests/test_examples_system.py

@@ -0,0 +1,25 @@
+import numpy as np
+
+from cryo_tank.config import P_WORKING, T_IN_HE
+from examples.cryo_tank_cylinder_system import run_system
+
+
+def test_cryo_tank_cylinder_system_passes_he_boundary_to_cylinder():
+    history = run_system(t_end=20.0, max_step=1.0)
+
+    assert len(history['t']) > 2
+    assert np.all(history['mdot_He'] >= 0.0)
+    assert np.allclose(history['mdot_tank_inlet'], history['mdot_He'])
+    assert np.allclose(history['mdot_cylinder_out'], history['mdot_He'])
+    assert np.allclose(history['T_tank_ullage'], history['T_ull'])
+    assert np.allclose(history['T_tank_liq'], history['T_liq'])
+    assert history['m_cylinder'][-1] < history['m_cylinder'][0]
+    assert history['U_cylinder'][-1] < history['U_cylinder'][0]
+    assert np.all(history['P_he_source'] > history['P_tank'])
+    assert np.allclose(history['P_he_source'], history['P_cylinder'])
+    assert np.allclose(history['P_he_tank_inlet'], history['P_tank'])
+    assert np.allclose(history['P_he_boundary'], P_WORKING)
+    assert np.allclose(history['T_he_tank_inlet'], T_IN_HE)
+    assert np.allclose(history['T_he_boundary'], T_IN_HE)
+    assert np.all(history['h_he_boundary'] > 0.0)
+    assert np.allclose(history['edot_He'], history['mdot_He'] * history['h_he_boundary'])

tests/test_gas_cylinder.py

@@ -0,0 +1,89 @@
+import pytest
+
+from cylinder.gas_cylinder import (
+    DEFAULT_FLUID,
+    DEFAULT_P_INIT,
+    DEFAULT_T_INIT,
+    DEFAULT_VOLUME,
+    HighPressureGasCylinder,
+)
+
+
+def test_default_high_pressure_cylinder_state_matches_inputs():
+    cylinder = HighPressureGasCylinder()
+
+    assert cylinder.fluid == DEFAULT_FLUID
+    assert cylinder.V == DEFAULT_VOLUME
+    assert abs(cylinder.P - DEFAULT_P_INIT) / DEFAULT_P_INIT < 1e-8
+    assert abs(cylinder.T - DEFAULT_T_INIT) < 1e-8
+    assert cylinder.mass > 0
+    assert cylinder.U > 0
+
+
+def test_cylinder_accepts_configurable_fluid_pressure_volume_temperature():
+    cylinder = HighPressureGasCylinder(
+        fluid="Nitrogen",
+        P_init=2.0e6,
+        V=50e-3,
+        T_init=300.0,
+    )
+
+    assert cylinder.fluid == "Nitrogen"
+    assert abs(cylinder.P - 2.0e6) / 2.0e6 < 1e-8
+    assert abs(cylinder.T - 300.0) < 1e-8
+    assert abs(cylinder.V - 50e-3) < 1e-15
+    assert cylinder.mass > 0
+
+
+def test_cylinder_state_summary_omits_unused_thermodynamic_helpers():
+    cylinder = HighPressureGasCylinder()
+    summary = cylinder.state_summary()
+
+    assert set(summary) == {
+        "fluid",
+        "V",
+        "mass",
+        "rho",
+        "P",
+        "T",
+        "U",
+        "u",
+        "h",
+    }
+
+
+def test_cylinder_ghost_state_is_stagnant_conservative_vector():
+    cylinder = HighPressureGasCylinder()
+    ghost = cylinder.ghost_state()
+
+    assert ghost.shape == (3,)
+    assert abs(ghost[0] - cylinder.rho) < 1e-12
+    assert ghost[1] == 0.0
+    assert abs(ghost[2] - cylinder.rho * cylinder.specific_internal_energy) < 1e-6
+
+
+def test_cylinder_apply_flux_outflow_reduces_mass_energy_and_pressure():
+    cylinder = HighPressureGasCylinder()
+    mass_before = cylinder.mass
+    U_before = cylinder.U
+    P_before = cylinder.P
+
+    mdot = 1e-3
+    edot = mdot * cylinder.h
+    dt = 1.0
+    cylinder.apply_flux(mdot=mdot, edot=edot, dt=dt, sign=-1)
+
+    assert abs(cylinder.mass - (mass_before - mdot * dt)) < 1e-12
+    assert abs(cylinder.U - (U_before - edot * dt)) < 1e-6
+    assert cylinder.P < P_before
+
+
+def test_cylinder_rejects_invalid_inputs():
+    with pytest.raises(ValueError):
+        HighPressureGasCylinder(fluid="")
+    with pytest.raises(ValueError):
+        HighPressureGasCylinder(P_init=0.0)
+    with pytest.raises(ValueError):
+        HighPressureGasCylinder(V=0.0)
+    with pytest.raises(ValueError):
+        HighPressureGasCylinder(T_init=0.0)

tests/test_integration.py

@@ -5,9 +5,9 @@ are the load-bearing tests for the "flux doubling" coupling mechanism.
 """
 import numpy as np
 import pytest
-from tank import Tank
+from tank_pipe.tank import Tank
-from pipe import Pipe
+from tank_pipe.pipe import Pipe
-from solver import run
+from tank_pipe.solver import run
 
 
 GAMMA = 1.4

tests/test_pipe.py

@@ -1,7 +1,7 @@
 # tests/test_pipe.py
 import numpy as np
 import pytest
-from pipe import Pipe
+from tank_pipe.pipe import Pipe
 
 
 GAMMA = 1.4

tests/test_riemann.py

@@ -1,6 +1,6 @@
 # tests/test_riemann.py
 import numpy as np
-from riemann import hll_flux
+from tank_pipe.riemann import hll_flux
 
 
 GAMMA = 1.4

tests/test_tank.py

@@ -1,7 +1,7 @@
 # tests/test_tank.py
 import numpy as np
 import pytest
-from tank import Tank
+from tank_pipe.tank import Tank
 
 
 GAMMA = 1.4