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