From f0310b6b404785210d8891a153f5f5bc5ed6d083 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=8D=A2=E4=BA=AC=E6=B3=BD?= Date: Wed, 10 Jun 2026 03:02:59 +0000 Subject: [PATCH] =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E8=81=94=E5=90=88=E4=BB=BF?= =?UTF-8?q?=E7=9C=9F=E5=86=85=E5=AE=B9=EF=BC=8C=E4=BF=AE=E6=94=B9index.md?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/__init__.py | 1 + examples/cryo_tank_cylinder_system.py | 389 ++++++++++++++++++++++++++ tests/test_examples_system.py | 25 ++ 3 files changed, 415 insertions(+) create mode 100644 examples/__init__.py create mode 100644 examples/cryo_tank_cylinder_system.py create mode 100644 tests/test_examples_system.py diff --git a/examples/__init__.py b/examples/__init__.py new file mode 100644 index 0000000..cf91e52 --- /dev/null +++ b/examples/__init__.py @@ -0,0 +1 @@ +"""Example system simulations.""" diff --git a/examples/cryo_tank_cylinder_system.py b/examples/cryo_tank_cylinder_system.py new file mode 100644 index 0000000..88d350a --- /dev/null +++ b/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() diff --git a/tests/test_examples_system.py b/tests/test_examples_system.py new file mode 100644 index 0000000..abf63de --- /dev/null +++ b/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'])