# HydDown hydrogen/other gas depressurisation
# Copyright (c) 2021-2025 Anders Andreasen
# Published under an MIT license
"""
Main calculation engine for pressure vessel filling and discharge simulations.
This module contains the HydDown class, which is the core of the HydDown package.
It integrates mass and energy balances over time to simulate pressure vessel
depressurization (discharge) and pressurization (filling) with heat transfer effects.
The HydDown class:
- Reads and validates YAML input defining vessel geometry, initial conditions,
calculation type, valve parameters, and heat transfer settings
- Initializes thermodynamic state using CoolProp for fluid properties
- Integrates mass and energy balances using explicit Euler time stepping
- Calculates heat transfer between fluid and vessel wall
- Handles various thermodynamic paths: isothermal, isenthalpic, isentropic,
constant internal energy, and full energy balance
- Supports multiple valve types: orifice, control valve, relief valve, constant mass flow
- Models fire heat loads using Stefan-Boltzmann approach
- Tracks two-phase systems with separate gas/liquid temperatures
- Can model 1-D transient heat conduction through vessel walls (composite materials)
- Stores time-series results and provides plotting capabilities
Calculation types:
- isothermal: Constant temperature (very slow process with large heat reservoir)
- isenthalpic: Constant enthalpy (adiabatic, no work)
- isentropic: Constant entropy (adiabatic with PV work)
- specified_U: Constant internal energy
- energybalance: Full energy balance with heat transfer and work
Heat transfer modes (for energybalance):
- fixed_U: Fixed overall heat transfer coefficient
- fixed_Q: Fixed heat input rate
- specified_h: Specified internal/external heat transfer coefficients
- detailed: 1-D transient conduction through vessel wall
- fire: External fire heat load using Stefan-Boltzmann equation
Valve types:
- orifice: Compressible flow through orifice (Yellow Book equation)
- control_valve: Control valve with Cv characteristic
- relief_valve: API 520/521 relief valve sizing
- mdot: Constant mass flow rate
The integration scheme uses explicit Euler method with user-specified time step.
Results are stored in numpy arrays for time, pressure, temperature, mass flow, etc.
Typical usage:
import yaml
from hyddown import HydDown
with open('input.yml') as f:
input_data = yaml.load(f, Loader=yaml.FullLoader)
hdown = HydDown(input_data)
hdown.run()
hdown.plot()
"""
import math
import numpy as np
import pandas as pd
from tqdm import tqdm
from scipy.optimize import minimize
from scipy.optimize import root_scalar
from CoolProp.CoolProp import PropsSI
import CoolProp.CoolProp as CP
from hyddown import transport as tp
from hyddown import validator
from hyddown import fire
from hyddown import thermesh as tm
import fluids
[docs]
class HydDown:
"""
Main class to to hold problem definition, running problem, storing results etc.
"""
[docs]
def __init__(self, input):
"""
Parameters
----------
input : dict
Dict holding problem definition
"""
self.input = input
self.verbose = 0
self.isrun = False
self.validate_input()
self.read_input()
self.initialize()
# sb_heat_load = np.loadtxt(
# "C:\\Users\\AndersAndreasen\\Documents\\GitHub\\HydDown\\src\\hyddown\\examples\\LPG-heat_load.txt"
# )
# self.sb_heat_load = lambda x: np.interp(
# x, sb_heat_load[:, 0], sb_heat_load[:, 1]
# )
[docs]
def initialize(self):
"""
Preparing for running problem by creating the fluid objects required
instantiating arrays for storing time-dependent results, setting additional
required class attributes.
"""
self.vol = self.inner_vol.V_total
self.vol_tot = self.outer_vol.V_total
self.vol_solid = self.vol_tot - self.vol
self.surf_area_outer = self.outer_vol.A
self.surf_area_inner = self.inner_vol.A
self.fluid = CP.AbstractState("HEOS", self.comp)
if "&" in self.comp:
self.fluid.specify_phase(CP.iphase_gas)
self.fluid.set_mole_fractions(self.molefracs)
self.transport_fluid = CP.AbstractState("HEOS", self.compSRK)
self.transport_fluid.specify_phase(CP.iphase_gas)
self.transport_fluid.set_mole_fractions(self.molefracs)
self.transport_fluid_wet = CP.AbstractState("HEOS", self.compSRK)
self.transport_fluid_wet.specify_phase(CP.iphase_liquid)
self.transport_fluid_wet.set_mole_fractions(self.molefracs)
self.vent_fluid = CP.AbstractState("HEOS", self.comp)
self.vent_fluid.specify_phase(CP.iphase_gas)
self.vent_fluid.set_mole_fractions(self.molefracs)
if "liquid_level" in self.input["vessel"]:
ll = self.input["vessel"]["liquid_level"]
V_liquid = self.inner_vol.V_from_h(ll)
self.fluid.update(CP.PQ_INPUTS, self.p0, 0)
liq_rho = self.fluid.rhomass()
m_liq = V_liquid * liq_rho
self.fluid.update(CP.PQ_INPUTS, self.p0, 1)
gas_rho = self.fluid.rhomass()
V_vapour = self.inner_vol.V_total - V_liquid
m_vap = V_vapour * gas_rho
m_tot = m_liq + m_vap
rho0 = m_tot / self.inner_vol.V_total
self.Q0 = m_vap / m_tot
self.fluid.update(CP.PQ_INPUTS, self.p0, self.Q0)
self.T0 = self.fluid.T()
self.liquid_level0 = ll
self.vent_fluid.update(CP.PQ_INPUTS, self.p0, 1.0)
else:
self.fluid.update(CP.PT_INPUTS, self.p0, self.T0)
self.liquid_level0 = 0.0
self.Q0 = self.fluid.Q()
self.vent_fluid.update(CP.PT_INPUTS, self.p0, self.T0)
self.res_fluid = CP.AbstractState("HEOS", self.comp)
self.res_fluid.set_mole_fractions(self.molefracs)
if self.input["valve"]["flow"] == "filling":
self.res_fluid.update(CP.PT_INPUTS, self.p_back, self.T0)
# Non-equilibrium model: Create separate AbstractState objects for gas and liquid
if self.non_equilibrium:
# Create parent phase objects - these CAN enter two-phase region
# Do NOT use specify_phase - we want to detect phase change via quality
self.fluid_gas = CP.AbstractState("HEOS", self.comp)
self.fluid_gas.set_mole_fractions(self.molefracs)
self.fluid_liquid = CP.AbstractState("HEOS", self.comp)
self.fluid_liquid.set_mole_fractions(self.molefracs)
# Initialize at saturation conditions (equilibrium at t=0)
if "liquid_level" in self.input["vessel"]:
# Two-phase initial condition
self.fluid_gas.update(CP.PQ_INPUTS, self.p0, 1.0)
self.fluid_liquid.update(CP.PQ_INPUTS, self.p0, 0.0)
self.T_gas0 = self.fluid_gas.T() # Saturated vapour
self.T_liquid0 = self.fluid_liquid.T() # Saturated liquid
# Calculate initial masses
ll = self.input["vessel"]["liquid_level"]
V_liquid = self.inner_vol.V_from_h(ll)
self.m_liquid0 = V_liquid * self.fluid_liquid.rhomass()
V_vapour = self.inner_vol.V_total - V_liquid
self.m_gas0 = V_vapour * self.fluid_gas.rhomass()
else:
# Single-phase gas initial condition
self.fluid_gas.update(CP.PT_INPUTS, self.p0, self.T0)
self.T_gas0 = self.T0
self.m_gas0 = self.fluid_gas.rhomass() * self.vol
self.m_liquid0 = 0.0
self.T_liquid0 = self.T0 # Not used if no liquid
# data storage
data_len = int(self.time_tot / self.tstep)
self.rho = np.zeros(data_len)
self.T_fluid = np.zeros(data_len)
self.T_vent = np.zeros(data_len)
self.T_vessel = np.zeros(data_len)
self.T_vessel_wetted = np.zeros(data_len)
self.T_inner_wall = np.zeros(data_len)
self.T_inner_wall_wetted = np.zeros(data_len)
self.T_outer_wall = np.zeros(data_len)
self.T_outer_wall_wetted = np.zeros(data_len)
self.T_bonded_wall = np.zeros(data_len)
self.T_bonded_wall_wetted = np.zeros(data_len)
self.Q_outer = np.zeros(data_len)
self.Q_inner = np.zeros(data_len)
self.Q_outer_wetted = np.zeros(data_len)
self.Q_inner_wetted = np.zeros(data_len)
self.q_outer = np.zeros(data_len)
self.q_inner = np.zeros(data_len)
self.q_outer_wetted = np.zeros(data_len)
self.q_inner_wetted = np.zeros(data_len)
self.h_inside = np.zeros(data_len)
self.h_inside_wetted = np.zeros(data_len)
self.h_gas_liquid = np.zeros(data_len) # Gas-liquid interfacial HTC for NEM
# Initialize Biot arrays as NaN (will be calculated if thermal_conductivity_biot specified)
self.Biot = np.full(
data_len, np.nan
) # Biot number for lumped capacitance validation
self.Biot_wetted = np.full(data_len, np.nan) # Biot number for wetted region
self.T_vent = np.zeros(data_len)
self.H_mass = np.zeros(data_len)
self.S_mass = np.zeros(data_len)
self.U_mass = np.zeros(data_len)
self.U_tot = np.zeros(data_len)
self.U_res = np.zeros(data_len)
self.P = np.zeros(data_len)
self.mass_fluid = np.zeros(data_len)
self.mass_rate = np.zeros(data_len)
self.time_array = np.zeros(data_len)
self.relief_area = np.zeros(data_len)
self.temp_profile = []
self.rho0 = self.fluid.rhomass()
self.m0 = self.rho0 * self.vol
self.MW = self.fluid.molar_mass()
self.vapour_mole_fraction = np.zeros(data_len)
self.vapour_mass_fraction = np.zeros(data_len)
self.vapour_volume_fraction = np.zeros(data_len)
self.liquid_level = np.zeros(data_len)
# Non-equilibrium model data storage
if self.non_equilibrium:
self.T_gas = np.zeros(data_len) # Gas phase temperature
self.T_liquid = np.zeros(data_len) # Liquid phase temperature
self.m_gas = np.zeros(data_len) # Gas phase mass
self.m_liquid = np.zeros(data_len) # Liquid phase mass
self.U_gas = np.zeros(data_len) # Gas phase internal energy (specific)
self.U_liquid = np.zeros(data_len) # Liquid phase internal energy (specific)
self.rho_gas = np.zeros(data_len) # Gas phase density
self.rho_liquid = np.zeros(data_len) # Liquid phase density
self.mdot_phase_transfer = np.zeros(data_len) # Mass transfer rate (condensation > 0, evaporation < 0)
[docs]
def calc_liquid_level(self):
"""
Calculate liquid level height based on current two-phase fluid state.
For two-phase systems (0 ≤ quality ≤ 1), calculates the height of liquid
phase in the vessel based on vapor quality, phase densities, and vessel geometry.
Uses vessel geometry from fluids.TANK to convert liquid volume to height.
Parameters
----------
fluid : CoolProp AbstractState
Current fluid state
Returns
-------
float
Liquid level height from vessel bottom [m].
Returns 0.0 for single-phase gas (quality > 1 or subcooled liquid).
"""
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
rho_liq = self.fluid.saturated_liquid_keyed_output(CP.iDmass)
rho_vap = self.fluid.saturated_vapor_keyed_output(CP.iDmass)
m_liq = self.fluid.rhomass() * self.inner_vol.V_total * (1 - self.fluid.Q())
V_liq = m_liq / rho_liq
h_liq = self.inner_vol.h_from_V(V_liq)
return h_liq
else:
return 0.0
[docs]
def PHres(self, T, P, H):
"""
Residual enthalpy function to be minimized during PH-problem.
Used by numerical optimizer (scipy.optimize.minimize) to find temperature
that satisfies constant pressure-enthalpy constraints for multicomponent
fluids. The optimizer adjusts T until residual approaches zero.
Parameters
----------
H : float
Enthalpy at initial/final conditions [J/kg]
P : float
Pressure at final conditions [Pa]
T : float
Updated estimate for the final temperature at (P,H) [K]
Returns
-------
float
Squared normalized enthalpy residual (dimensionless).
Zero when correct temperature is found.
"""
# Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars)
T_scalar = float(T.item()) if hasattr(T, "item") else float(T)
self.vent_fluid.update(CP.PT_INPUTS, P, T_scalar)
return ((H - self.vent_fluid.hmass()) / H) ** 2
[docs]
def PHres_relief(self, T, P, H):
"""
Residual enthalpy function for PH-problem during relief valve calculations.
Used by numerical optimizer (scipy.optimize.root_scalar) to find temperature
for relief valve discharge calculations. Similar to PHres() but uses the
main fluid state instead of vent_fluid state.
Parameters
----------
H : float
Enthalpy at initial/final conditions [J/kg]
P : float
Pressure at final conditions [Pa]
T : float
Updated estimate for the final temperature at (P,H) [K]
Returns
-------
float
Normalized enthalpy residual (dimensionless).
Zero when correct temperature is found.
"""
# Extract scalar from array (scipy optimizers pass arrays, CoolProp needs scalars)
T_scalar = float(T.item()) if hasattr(T, "item") else float(T)
self.fluid.update(CP.PT_INPUTS, P, T_scalar)
return (H - self.fluid.hmass()) / H
[docs]
def PHproblem(self, H, P, Tguess, relief=False):
"""
Defining a constant pressure, constant enthalpy problem i.e. typical adiabatic
problem like e.g. valve flow for the vented flow (during discharge).
For multicomponent mixture the final temperature is changed/optimised until the residual
enthalpy is near zero in an optimisation step. For single component fluid the coolprop
built in methods are used for speed.
Parameters
----------
H : float
Enthalpy at initial/final conditions
P : float
Pressure at final conditions.
Tguess : float
Initial guess for the final temperature at P,H
"""
# Multicomponent case
if "&" in self.species:
x0 = Tguess
if relief == False:
res = minimize(
self.PHres,
x0,
args=(P, H),
method="Nelder-Mead",
options={"xatol": 0.1, "fatol": 0.001},
)
T1 = res.x[0]
else:
res = root_scalar(
self.PHres_relief,
args=(P, H),
x0=x0,
method="newton",
)
T1 = res.root
# single component fluid case
else:
T1 = PropsSI("T", "P", P, "H", H, self.species)
return T1
[docs]
def UDres(self, x, U, rho):
"""
Residual U-rho to be minimised during a U-rho/UV-problem
Parameters
----------
U : float
Internal energy at final conditions
rho : float
Density at final conditions
"""
self.fluid.update(CP.PT_INPUTS, x[0], x[1])
return ((U - self.fluid.umass()) / U) ** 2 + (
(rho - self.fluid.rhomass()) / rho
) ** 2
[docs]
def UDproblem(self, U, rho, Pguess, Tguess):
"""
Defining a constant UV problem i.e. constant internal energy and density/volume
problem relevant for the 1. law of thermodynamics.
For multicomponent mixture the final temperature/pressure is changed/optimised until the
residual U/rho is near zero. For single component fluid the coolprop
built in methods are used for speed.
Parameters
----------
U : float
Internal energy at final conditions
rho : float
Density at final conditions.
Pguess : float
Initial guess for the final pressure at U, rho
Tguess : float
Initial guess for the final temperature at U, rho
"""
if "&" in self.species:
x0 = [Pguess, Tguess]
res = minimize(
self.UDres,
x0,
args=(U, rho),
method="Nelder-Mead",
options={"xatol": 0.1, "fatol": 0.001},
)
P1 = res.x[0]
T1 = res.x[1]
Ures = U - self.fluid.umass()
else:
P1 = PropsSI("P", "D", rho, "U", U, self.species)
T1 = PropsSI("T", "D", rho, "U", U, self.species)
Ures = 0
return P1, T1, Ures
[docs]
def nem_objective_function(self, x, m_gas, m_liquid, U_gas_tot, U_liquid_tot, T_gas_prev, T_liquid_prev):
"""
Objective function for NEM solver (inspired by ORS-openthermo).
Solves for (T_gas, T_liquid, P) simultaneously by minimizing
normalized residuals for volume and energy constraints.
Parameters
----------
x : array
[T_gas, T_liquid, P]
m_gas : float
Gas mass [kg]
m_liquid : float
Liquid mass [kg]
U_gas_tot : float
Total gas internal energy [J]
U_liquid_tot : float
Total liquid internal energy [J]
T_gas_prev : float
Previous gas temperature (for bounds checking) [K]
T_liquid_prev : float
Previous liquid temperature (for bounds checking) [K]
Returns
-------
float
Objective value (sum of normalized squared residuals)
"""
T_gas, T_liquid, P = x
# Sanity checks
if P <= 0 or T_gas <= 0 or T_liquid <= 0:
return 1e10
# Sanity checks on inputs - return high penalty if out of reasonable range
if P <= 1e4 or P > 1e8: # 0.1 to 1000 bar
return 1e10
if T_gas <= 100 or T_gas > 1000: # 100 K to 1000 K
return 1e10
if T_liquid <= 100 or T_liquid > 1000:
return 1e10
# Update gas state with T, P
# Handle case where PT is at saturation (CoolProp will fail)
if m_gas > 1e-6:
try:
self.fluid_gas.update(CP.PT_INPUTS, P, T_gas)
rho_gas = self.fluid_gas.rhomass()
U_gas_calc = self.fluid_gas.umass() * m_gas
V_gas = m_gas / rho_gas
except ValueError as e:
# PT failed - likely at or near saturation
# Penalize this state heavily to push optimizer away from saturation line
return 1e10
else:
U_gas_calc = 0.0
V_gas = 0.0
# Update liquid state with T, P
if m_liquid > 1e-6:
try:
self.fluid_liquid.update(CP.PT_INPUTS, P, T_liquid)
rho_liquid = self.fluid_liquid.rhomass()
U_liquid_calc = self.fluid_liquid.umass() * m_liquid
V_liquid = m_liquid / rho_liquid
except ValueError as e:
# PT failed - likely at or near saturation
# Penalize this state heavily
return 1e10
else:
U_liquid_calc = 0.0
V_liquid = 0.0
# Calculate normalized squared residuals
V_total = V_gas + V_liquid
# Volume constraint (normalized by vessel volume)
vol_res = ((V_total - self.vol) / self.vol) ** 2
# Energy constraints (normalized by target energy)
if m_gas > 1e-6 and U_gas_tot > 0:
energy_res_gas = ((U_gas_calc - U_gas_tot) / U_gas_tot) ** 2
else:
energy_res_gas = 0.0
if m_liquid > 1e-6 and U_liquid_tot > 0:
energy_res_liq = ((U_liquid_calc - U_liquid_tot) / U_liquid_tot) ** 2
else:
energy_res_liq = 0.0
# Combined objective with equal weights
# All residuals are normalized squared differences
objective = vol_res + energy_res_gas + energy_res_liq
# Debug: Store individual residuals for diagnosis (only if requested)
if hasattr(self, '_debug_nem_residuals') and self._debug_nem_residuals:
print(f" P={P/1e5:.2f} bar, T_g={T_gas:.1f} K, T_l={T_liquid:.1f} K")
print(f" vol_res={vol_res:.2e}, U_gas_res={energy_res_gas:.2e}, U_liq_res={energy_res_liq:.2e}")
print(f" V_calc={V_total:.4f} m³ vs V_target={self.vol:.4f} m³")
print(f" U_gas_calc={U_gas_calc/1e6:.3f} MJ vs U_gas_target={U_gas_tot/1e6:.3f} MJ")
print(f" U_liq_calc={U_liquid_calc/1e6:.3f} MJ vs U_liq_target={U_liquid_tot/1e6:.3f} MJ")
return objective
[docs]
def nem_residual(self, P, m_gas, m_liquid, U_gas_tot, U_liquid_tot):
"""
Residual function for non-equilibrium model (NEM) solver.
Solves for common pressure P such that total volume equals vessel volume.
Key insight: Both phases are in mechanical equilibrium (same pressure P).
Given P and U for each phase, CoolProp gives us rho and T via DmassUmass inputs.
The constraint is: V_gas(P,U_gas) + V_liquid(P,U_liquid) = V_vessel
Parameters
----------
P : float
Common pressure for both phases [Pa]
m_gas : float
Gas phase mass [kg]
m_liquid : float
Liquid phase mass [kg]
U_gas_tot : float
Total gas phase internal energy [J]
U_liquid_tot : float
Total liquid phase internal energy [J]
Returns
-------
float
Volume residual (V_total - V_vessel) / V_vessel
"""
# Ensure positive pressure
if P <= 0:
return 1e10
V_total = 0.0
# Gas phase: given P and U, find rho (and T implicitly)
# We need to iterate to find rho_gas such that fluid_gas(rho_gas, U_gas) gives pressure P
if m_gas > 1e-6:
U_gas_spec = U_gas_tot / m_gas
try:
# For given P and U, find the density
# This requires iteration since we can't directly specify P,U to CoolProp
# We use PUmass_INPUTS if available, otherwise iterate
self.fluid_gas.update(CP.PUmass_INPUTS, P, U_gas_spec)
rho_gas = self.fluid_gas.rhomass()
V_gas = m_gas / rho_gas
V_total += V_gas
except:
# PUmass might not work for all fluids, try iteration
# Initial guess: use ideal gas for first estimate
from scipy.optimize import brentq
def rho_residual(rho):
try:
self.fluid_gas.update(CP.DmassUmass_INPUTS, rho, U_gas_spec)
return self.fluid_gas.p() - P
except:
return 1e10
try:
# Search for density that gives the target pressure
rho_gas = brentq(rho_residual, 0.01, 1000.0, xtol=1e-6)
V_gas = m_gas / rho_gas
V_total += V_gas
except:
return 1e10
# Liquid phase: given P and U, find rho (and T implicitly)
if m_liquid > 1e-6:
U_liquid_spec = U_liquid_tot / m_liquid
try:
self.fluid_liquid.update(CP.PUmass_INPUTS, P, U_liquid_spec)
rho_liquid = self.fluid_liquid.rhomass()
V_liquid = m_liquid / rho_liquid
V_total += V_liquid
except:
# PUmass might not work, try iteration
from scipy.optimize import brentq
def rho_residual(rho):
try:
self.fluid_liquid.update(CP.DmassUmass_INPUTS, rho, U_liquid_spec)
return self.fluid_liquid.p() - P
except:
return 1e10
try:
# Search for density that gives the target pressure
rho_liquid = brentq(rho_residual, 100.0, 2000.0, xtol=1e-6)
V_liquid = m_liquid / rho_liquid
V_total += V_liquid
except:
return 1e10
# Return volume residual
vol_residual = (V_total - self.vol) / self.vol
return vol_residual
[docs]
def nem_solve_state(self, m_gas, m_liquid, U_gas_tot, U_liquid_tot, P_guess, T_gas_guess, T_liquid_guess):
"""
Solve for (T_gas, T_liquid, P) in non-equilibrium model.
Uses optimization approach inspired by ORS-openthermo: minimize combined
residuals for volume and energy constraints.
Parameters
----------
m_gas : float
Gas phase mass [kg]
m_liquid : float
Liquid phase mass [kg]
U_gas_tot : float
Total gas internal energy [J]
U_liquid_tot : float
Total liquid internal energy [J]
P_guess : float
Initial pressure guess [Pa]
T_gas_guess : float
Initial gas temperature guess [K]
T_liquid_guess : float
Initial liquid temperature guess [K]
Returns
-------
tuple
(P, T_gas, T_liquid) - solved state [Pa, K, K]
"""
from scipy.optimize import minimize
# Initial guess
x0 = [T_gas_guess, T_liquid_guess, P_guess]
# Bounds: allow reasonable deviations from previous state
bounds = [
(max(T_gas_guess - 50, 200), T_gas_guess + 100), # T_gas bounds
(max(T_liquid_guess - 50, 200), T_liquid_guess + 100), # T_liquid bounds
(max(P_guess * 0.5, 1e5), P_guess * 2.0) # P bounds
]
# Use Nelder-Mead (doesn't require derivatives, works well for non-smooth problems)
result = minimize(
fun=self.nem_objective_function,
x0=x0,
args=(m_gas, m_liquid, U_gas_tot, U_liquid_tot, T_gas_guess, T_liquid_guess),
method='Nelder-Mead',
bounds=bounds,
options={'maxiter': 1000, 'xatol': 1e-6, 'fatol': 1e-9}
)
# Accept solution if objective is reasonable
# Note: Objective = vol_res² + U_gas_res² + U_liq_res² (all normalized)
# Objective < 1e-4 means ~0.01% error in each constraint
# Objective < 1e-2 means ~1% error in each constraint
if result.fun < 1e-2:
T_gas, T_liquid, P = result.x
return P, T_gas, T_liquid
# If convergence is poor, print warning but still return result
# This can happen when constraints are difficult to satisfy exactly
if result.fun < 0.1: # ~10% error - still usable
print(f"Warning: NEM solver converged with objective={result.fun:.2e} at P={P_guess/1e5:.1f} bar")
T_gas, T_liquid, P = result.x
return P, T_gas, T_liquid
# If Nelder-Mead didn't converge well enough, raise error
raise RuntimeError(f"NEM solver failed. Objective={result.fun:.2e}, P_guess={P_guess/1e5:.2f} bar, T_gas={T_gas_guess:.1f} K, T_liq={T_liquid_guess:.1f} K")
[docs]
def nem_calc_phase_transfer(self, P, T_gas, T_liquid, m_gas, m_liquid, dt):
"""
Calculate phase transfer based on vapor quality from CoolProp thermodynamic updates.
After updating each phase with its internal energy, CoolProp returns a vapor quality:
- For liquid phase: quality > 0 means some liquid has evaporated (transfer to gas)
- For gas phase: quality < 1 means some gas has condensed (transfer to liquid)
This approach uses thermodynamics directly rather than relaxation assumptions.
Parameters
----------
P : float
Current pressure [Pa]
T_gas : float
Gas phase temperature [K]
T_liquid : float
Liquid phase temperature [K]
m_gas : float
Gas phase mass [kg]
m_liquid : float
Liquid phase mass [kg]
dt : float
Time step [s]
Returns
-------
dm_transfer : float
Net mass transfer [kg] (positive = condensation gas→liquid, negative = evaporation liquid→gas)
Note: This function signature kept for compatibility, but actual phase transfer
is now calculated directly in the main integration loop using vapor quality.
"""
# This function is now a placeholder - actual phase transfer calculated
# in main loop using vapor quality from CoolProp updates
return 0.0
[docs]
def run(self, disable_pbar=True):
"""
Routine for running the actual problem defined i.e. integrating the mass and energy balances
"""
# Inititialise / setting initial values for t=0
if self.isrun is True:
self.initialize()
input = self.input
self.rho[0] = self.rho0
self.T_fluid[0] = self.T0
self.T_vessel[0] = self.T0
self.T_inner_wall[0] = self.T0
self.T_outer_wall[0] = self.T0
self.T_vessel_wetted[0] = self.T0
self.T_inner_wall_wetted[0] = self.T0
self.T_outer_wall_wetted[0] = self.T0
self.T_bonded_wall[0] = self.T0
self.T_bonded_wall_wetted[0] = self.T0
self.liquid_level[0] = self.liquid_level0
if self.input["valve"]["flow"] == "discharge":
self.T_vent[0] = self.T0
self.H_mass[0] = self.fluid.hmass()
self.S_mass[0] = self.fluid.smass()
self.U_mass[0] = self.fluid.umass()
self.U_tot[0] = self.fluid.umass() * self.m0
self.P[0] = self.p0
self.mass_fluid[0] = self.m0
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
self.vapour_mole_fraction[0] = self.fluid.Q()
else:
self.vapour_mole_fraction[0] = 1.0
# Non-equilibrium model initialization
if self.non_equilibrium:
self.T_gas[0] = self.T_gas0
self.T_liquid[0] = self.T_liquid0
self.m_gas[0] = self.m_gas0
self.m_liquid[0] = self.m_liquid0
if self.m_liquid0 > 0:
# Two-phase: try saturation, fall back to +1°C superheat
self.fluid_liquid.update(CP.PQ_INPUTS, self.p0, 0.0)
self.U_liquid[0] = self.fluid_liquid.umass()
self.rho_liquid[0] = self.fluid_liquid.rhomass()
try:
self.fluid_gas.update(CP.PQ_INPUTS, self.p0, 1.0)
_rho = self.fluid_gas.rhomass()
_U = self.fluid_gas.umass()
# Verify DmassUmass round-trip at saturation
self.fluid_gas.update(CP.DmassUmass_INPUTS, _rho, _U)
# Verify quality is valid (Q=1 at saturation boundary)
_Q = self.fluid_gas.Q()
if _Q < 0 or _Q > 1:
raise ValueError("Invalid quality at saturation")
except Exception:
# Saturation boundary unstable - add slight superheat
self.T_gas0 += 1.0
self.T_gas[0] = self.T_gas0
self.fluid_gas.update(CP.PT_INPUTS, self.p0, self.T_gas0)
else:
# Single-phase gas
self.fluid_gas.update(CP.PT_INPUTS, self.p0, self.T_gas0)
self.U_liquid[0] = 0.0
self.rho_liquid[0] = 0.0
self.U_gas[0] = self.fluid_gas.umass()
self.rho_gas[0] = self.fluid_gas.rhomass()
self.mdot_phase_transfer[0] = 0.0
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
cpcv = self.fluid.saturated_vapor_keyed_output(CP.iCpmolar) / (
self.fluid.saturated_vapor_keyed_output(CP.iCpmolar) - 8.314
)
rho0 = self.fluid.saturated_vapor_keyed_output(CP.iDmass)
Z = self.fluid.saturated_vapor_keyed_output(CP.iZ)
else:
cpcv = self.fluid.cp0molar() / (self.fluid.cp0molar() - 8.314)
rho0 = self.rho0
Z = self.fluid.compressibility_factor()
massflow_stop_switch = 0
# Calculating initial mass rate for t=0 depending on mass flow device
# and filling/discharge mode
if input["valve"]["type"] == "orifice":
if input["valve"]["flow"] == "filling":
k = self.res_fluid.cp0molar() / (self.res_fluid.cp0molar() - 8.314)
self.mass_rate[0] = -tp.gas_release_rate(
self.p_back,
self.p0,
self.res_fluid.rhomass(),
k,
self.CD,
self.D_orifice**2 / 4 * math.pi,
)
else:
self.mass_rate[0] = tp.gas_release_rate(
self.p0,
self.p_back,
rho0,
cpcv,
self.CD,
self.D_orifice**2 / 4 * math.pi,
)
elif input["valve"]["type"] == "hem_release":
if input["valve"]["flow"] == "filling":
raise ValueError(
"Unsupported valve: ",
input["valve"]["type"],
" for vessel filling.",
)
else:
self.mass_rate[0] = tp.hem_release_rate(
self.p0,
self.p_back,
self.CD,
self.D_orifice**2 / 4 * math.pi,
self.fluid,
)
elif input["valve"]["type"] == "mdot":
if "mdot" in input["valve"].keys() and "time" in input["valve"].keys():
mdot = np.asarray(input["valve"]["mdot"])
time = np.asarray(input["valve"]["time"])
max_i = int(time[-1] / self.tstep)
interp_time = np.linspace(
0,
self.tstep * len(self.time_array),
len(self.time_array),
endpoint=False,
)[:max_i]
self.mass_rate[:max_i] = np.interp(interp_time, time, mdot)
if input["valve"]["flow"] == "filling":
self.mass_rate *= -1
else:
if input["valve"]["flow"] == "filling":
self.mass_rate[:] = -input["valve"]["mdot"]
else:
self.mass_rate[:] = input["valve"]["mdot"]
elif input["valve"]["type"] == "controlvalve":
Cv = tp.cv_vs_time(
self.Cv, 0, self.valve_time_constant, self.valve_characteristic
)
if input["valve"]["flow"] == "filling":
Z = self.res_fluid.compressibility_factor()
MW = self.MW
k = self.res_fluid.cp0molar() / (self.res_fluid.cp0molar() - 8.314)
self.mass_rate[0] = -tp.control_valve(
self.p_back, self.p0, self.T0, Z, MW, k, Cv
)
else:
MW = self.MW
k = cpcv
self.mass_rate[0] = tp.control_valve(
self.p0, self.p_back, self.T0, Z, MW, k, Cv
)
elif input["valve"]["type"] == "psv":
if input["valve"]["flow"] == "filling":
raise ValueError(
"Unsupported valve: ",
input["valve"]["type"],
" for vessel filling.",
)
self.mass_rate[0] = tp.relief_valve(
self.p0,
self.p_back,
self.Pset,
self.blowdown,
cpcv,
self.CD,
self.T0,
Z,
self.MW,
self.D_orifice**2 / 4 * math.pi,
)
self.time_array[0] = 0
# ============================================================================
# MAIN TIME INTEGRATION LOOP
# ============================================================================
# Integrate mass and energy balances forward in time using explicit Euler method
# At each time step:
# 1. Update mass inventory from mass flow rate
# 2. Calculate new density (mass/volume)
# 3. Determine thermodynamic state based on calculation method:
# - Simple methods (isenthalpic/isentropic/isothermal/constantU):
# Direct CoolProp update with (density, H/S/T/U)
# - Energy balance method: Solve for temperature from energy balance
# accounting for heat transfer, work, and enthalpy changes
# 4. Calculate heat transfer (if energybalance method)
# 5. Update vessel wall temperatures
# 6. Calculate mass flow rate for next time step
# ============================================================================
# Initialize heat transfer variables
T_profile, T_profile2 = 0, 0 # Temperature profiles for detailed wall model
relief_area = [] # Track relief valve area changes
for i in tqdm(
range(1, len(self.time_array)),
desc="hyddown",
disable=disable_pbar,
total=len(self.time_array),
):
self.time_array[i] = self.time_array[i - 1] + self.tstep
self.mass_fluid[i] = (
self.mass_fluid[i - 1] - self.mass_rate[i - 1] * self.tstep
)
self.rho[i] = self.mass_fluid[i] / self.vol
# ------------------------------------------------------------------------
# THERMODYNAMIC STATE UPDATE
# ------------------------------------------------------------------------
# For simple methods, density changes but one other property remains constant
# CoolProp can directly calculate (T, P) from (density, H/S/U/T) for single components
# For multicomponent fluids, numerical optimization is required (slower)
if self.method == "isenthalpic":
self.fluid.update(CP.DmassHmass_INPUTS, self.rho[i], self.H_mass[i - 1])
self.T_fluid[i] = self.fluid.T()
self.P[i] = self.fluid.p()
elif self.method == "isentropic":
self.fluid.update(CP.DmassSmass_INPUTS, self.rho[i], self.S_mass[i - 1])
self.T_fluid[i] = self.fluid.T()
self.P[i] = self.fluid.p()
elif self.method == "isothermal":
self.fluid.update(CP.DmassT_INPUTS, self.rho[i], self.T0)
self.T_fluid[i] = self.T0
self.P[i] = self.fluid.p()
elif self.method == "constantU":
self.fluid.update(CP.DmassUmass_INPUTS, self.rho[i], self.U_mass[i - 1])
self.T_fluid[i] = self.fluid.T()
self.P[i] = self.fluid.p()
# ------------------------------------------------------------------------
# ENERGY BALANCE METHOD (Most Complex Case)
# ------------------------------------------------------------------------
# Full energy balance accounting for:
# - Heat transfer between fluid and vessel wall (convection)
# - Heat transfer between vessel and environment (convection, radiation, fire)
# - Enthalpy change from mass flow (inlet/outlet streams)
# - Work done by expanding/compressing fluid
# - 1-D transient conduction through vessel wall (if detailed method)
#
# Energy balance on fluid:
# dU/dt = Q_transfer + H_inlet*dm_inlet/dt - H_outlet*dm_outlet/dt
#
# where U = internal energy, Q = heat transfer rate, H = specific enthalpy
elif self.method == "energybalance":
# ====================================================================
# HEAT TRANSFER COEFFICIENT CALCULATIONS
# ====================================================================
# Calculate convective heat transfer coefficients for specified_h or detailed methods
if self.heat_method == "specified_h" or self.heat_method == "detailed" or self.heat_method == "specified_q":
if self.h_in == "calc":
if self.vessel_orientation == "horizontal":
L = self.diameter
else:
L = self.length
if input["valve"]["flow"] == "filling":
# T_film = (self.T_fluid[i - 1] + self.T_vessel[i - 1]) / 2
# For NEM: use gas temperature for unwetted (gas-side) heat transfer
# For equilibrium: use bulk fluid temperature
if self.non_equilibrium and hasattr(self, 'T_gas'):
T_for_gas_side_htc = self.T_gas[i - 1]
else:
T_for_gas_side_htc = self.T_fluid[i - 1]
T_film = (
T_for_gas_side_htc + self.T_inner_wall[i - 1]
) / 2
self.transport_fluid.update(
CP.PT_INPUTS, self.P[i - 1], T_film
)
hi = tp.h_inside_mixed(
L,
# self.T_vessel[i - 1],
self.T_inner_wall[i - i],
T_for_gas_side_htc,
self.transport_fluid,
self.mass_rate[i - 1],
self.diameter,
)
# NEM: Use liquid properties for wetted heat transfer
# Equilibrium: Use existing logic with equilibrium fluid
if self.non_equilibrium:
# For NEM, if liquid exists, use nucleate boiling correlation
# No need to check quality - liquid phase is always liquid
liquid_exists = self.m_liquid[i-1] > 1e-6
if liquid_exists:
# Use liquid temperature for film temperature
T_film_wet = (self.T_liquid[i-1] + self.T_inner_wall_wetted[i-1]) / 2
try:
self.transport_fluid_wet.update(CP.PT_INPUTS, self.P[i-1], T_film_wet)
except:
# If film temperature fails, use liquid temperature
try:
self.transport_fluid_wet.update(CP.PT_INPUTS, self.P[i-1], self.T_liquid[i-1])
except:
# Last resort: use saturation
self.transport_fluid_wet.update(CP.PQ_INPUTS, self.P[i-1], 0.0)
# Update fluid_liquid to saturated state for h_inside_wetted
# (needed for surface tension and saturated properties)
self.fluid_liquid.update(CP.PQ_INPUTS, self.P[i-1], 0.0)
hiw = tp.h_inside_wetted(
L,
self.T_inner_wall_wetted[i - 1],
self.T_liquid[i - 1], # Use liquid temperature
self.transport_fluid_wet,
self.fluid_liquid, # Use liquid phase object at saturation
)
else:
# No liquid - use gas-side coefficient
hiw = hi
else:
# Equilibrium mode: use existing logic
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
self.transport_fluid_wet.update(
CP.PT_INPUTS, self.P[i - 1], T_film
)
hiw = tp.h_inside_wetted(
L,
self.T_inner_wall_wetted[i - 1],
self.T_fluid[i - 1],
self.transport_fluid_wet,
self.fluid,
)
else:
hiw = hi
else:
# For NEM: use gas temperature for unwetted (gas-side) heat transfer
# For equilibrium: use bulk fluid temperature
if self.non_equilibrium and hasattr(self, 'T_gas'):
T_for_gas_side_htc = self.T_gas[i - 1]
else:
T_for_gas_side_htc = self.T_fluid[i - 1]
T_film = (
T_for_gas_side_htc + self.T_inner_wall[i - 1]
) / 2
try:
self.transport_fluid.update(
CP.PT_INPUTS, self.P[i - 1], T_film
)
# Update transport fluid for wetted surface
# For NEM, will be updated later with liquid temperature
if not self.non_equilibrium:
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
self.transport_fluid_wet.update(
CP.PT_INPUTS, self.P[i - 1], self.T_fluid[i - 1]
)
except:
self.transport_fluid.update(
CP.PQ_INPUTS, self.P[i - 1], 1.0
)
hi = tp.h_inside(
L,
self.T_inner_wall[i - 1],
T_for_gas_side_htc,
self.transport_fluid,
)
# NEM: Check liquid phase for boiling, use liquid properties
# Equilibrium: Use existing logic with equilibrium fluid
if self.non_equilibrium:
# For NEM, if liquid exists, use nucleate boiling correlation
# No need to check quality - liquid phase is always liquid
liquid_exists = self.m_liquid[i-1] > 1e-6
if liquid_exists:
# Use liquid temperature for film temperature
T_film_wet = (self.T_liquid[i-1] + self.T_inner_wall_wetted[i-1]) / 2
try:
self.transport_fluid_wet.update(CP.PT_INPUTS, self.P[i-1], T_film_wet)
except:
# If film temperature fails, use liquid temperature
try:
self.transport_fluid_wet.update(CP.PT_INPUTS, self.P[i-1], self.T_liquid[i-1])
except:
# Last resort: use saturation
self.transport_fluid_wet.update(CP.PQ_INPUTS, self.P[i-1], 0.0)
# Update fluid_liquid to saturated state for h_inside_wetted
# (needed for surface tension and saturated properties)
self.fluid_liquid.update(CP.PQ_INPUTS, self.P[i-1], 0.0)
hiw = tp.h_inside_wetted(
L,
self.T_inner_wall_wetted[i - 1],
self.T_liquid[i - 1], # Use liquid temperature
self.transport_fluid_wet,
self.fluid_liquid, # Use liquid phase object at saturation
)
else:
# No liquid - use gas-side coefficient
hiw = hi
else:
# Equilibrium mode: use existing logic
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
hiw = tp.h_inside_wetted(
L,
self.T_inner_wall_wetted[i - 1],
self.T_fluid[i - 1],
self.transport_fluid_wet,
self.fluid,
)
else:
hiw = hi
else:
hi = self.h_in
hiw = self.h_in # Use same coefficient for wetted surface
self.h_inside[i] = hi
self.h_inside_wetted[i] = hiw
# ================================================================
# TWO-PHASE HEAT TRANSFER (Wetted vs Unwetted Areas)
# ================================================================
# For two-phase systems, split vessel into two regions:
# 1. Wetted area: Wall in contact with liquid phase
# 2. Unwetted area: Wall in contact with gas/vapor phase
# Different heat transfer coefficients and wall temperatures
# are calculated for each region based on liquid level
wetted_area = self.inner_vol.SA_from_h(self.liquid_level[i - 1])
if np.isnan(wetted_area):
wetted_area = 0
# Calculate outer wetted area for correct heat transfer on outer surface
# Add wall thickness to liquid level height since outer vessel reference
# is at bottom of wall, not bottom of inner surface
liquid_level_outer = self.liquid_level[i - 1] + self.thickness
wetted_area_outer = self.outer_vol.SA_from_h(liquid_level_outer)
if np.isnan(wetted_area_outer):
wetted_area_outer = 0
# Heat transfer from unwetted wall (gas side) to fluid
# For NEM: use gas temperature
# For equilibrium: use bulk fluid temperature
if self.non_equilibrium and hasattr(self, 'T_gas'):
T_for_gas_side_htc = self.T_gas[i - 1]
else:
T_for_gas_side_htc = self.T_fluid[i - 1]
self.Q_inner[i] = (
(self.surf_area_inner - wetted_area)
* hi
* (self.T_inner_wall[i - 1] - T_for_gas_side_htc)
)
self.q_inner[i] = hi * (
self.T_inner_wall[i - 1] - T_for_gas_side_htc
)
# Heat transfer from wetted wall (liquid side) to fluid
# Uses different heat transfer coefficient (hiw) for liquid contact
# For NEM: use liquid temperature
# For equilibrium: use bulk fluid temperature
if self.non_equilibrium and hasattr(self, 'T_liquid'):
T_fluid_wet = self.T_liquid[i - 1]
else:
T_fluid_wet = self.T_fluid[i - 1]
self.Q_inner_wetted[i] = (
wetted_area
* hiw
* (self.T_inner_wall_wetted[i - 1] - T_fluid_wet)
)
# Avoid division by zero when fully vapor (wetted_area = 0)
if wetted_area > 0:
self.q_inner_wetted[i] = self.Q_inner_wetted[i] / wetted_area
else:
self.q_inner_wetted[i] = 0.0
if np.isnan(self.Q_inner_wetted[i]):
self.Q_inner_wetted[i] = 0
# Heat transfer from environment to unwetted outer wall
# Use outer surface area directly (outer surface is exposed to environment)
if self.heat_method == "specified_q":
# User-defined heat flux (time-dependent or constant)
q_ext = self.q_outer_func(self.time_array[i])
self.Q_outer[i] = q_ext * (self.surf_area_outer - wetted_area_outer)
self.q_outer[i] = q_ext
self.Q_outer_wetted[i] = q_ext * wetted_area_outer
self.q_outer_wetted[i] = q_ext
else:
# specified_h or detailed: convective heat transfer
self.Q_outer[i] = (
(self.surf_area_outer - wetted_area_outer)
* self.h_out
* (self.Tamb - self.T_outer_wall[i - 1])
)
self.q_outer[i] = self.h_out * (
self.Tamb - self.T_outer_wall[i - 1]
)
self.Q_outer_wetted[i] = (
wetted_area_outer
* self.h_out
* (self.Tamb - self.T_outer_wall_wetted[i - 1])
)
self.q_outer_wetted[i] = self.h_out * (
self.Tamb - self.T_outer_wall_wetted[i - 1]
)
if np.isnan(self.Q_outer_wetted[i]):
self.Q_outer_wetted[i] = 0
self.T_vessel[i] = self.T_vessel[i - 1] + (
self.Q_outer[i] - self.Q_inner[i]
) * self.tstep / (
self.vessel_cp
* self.vessel_density
* self.vol_solid
* (self.inner_vol.A - wetted_area)
/ self.inner_vol.A
)
# Update wetted wall temperature only if wetted area exists
if wetted_area < 1e-10 or (self.Q_outer_wetted[i] == 0 and self.Q_inner_wetted[i] == 0):
self.T_vessel_wetted[i] = self.T_vessel_wetted[i - 1]
else:
self.T_vessel_wetted[i] = self.T_vessel_wetted[i - 1] + (
self.Q_outer_wetted[i] - self.Q_inner_wetted[i]
) * self.tstep / (
self.vessel_cp
* self.vessel_density
* self.vol_solid
* wetted_area
/ self.inner_vol.A
)
# ================================================================
# 1-D TRANSIENT HEAT CONDUCTION (Detailed Thermal Model)
# ================================================================
# Solve transient heat conduction through vessel wall using
# finite element method (thermesh module).
# Used for vessels with low thermal conductivity (Type III/IV)
# or composite materials where temperature gradient through
# wall thickness is significant.
#
# For single-layer walls: Uses isothermal material model
# For multi-layer walls: Uses piecewise linear model for liner+shell
#
# Time integration: Crank-Nicolson (theta=0.5) for stability
# Spatial discretization: Linear finite elements with 11 nodes
if "thermal_conductivity" in self.input["vessel"].keys():
theta = 0.5 # Crank-Nicolson scheme (unconditionally stable, 2nd order)
dt = (
self.tstep / 10
) # Sub-step for thermal solver (finer time resolution)
k, rho, cp = (
self.input["vessel"]["thermal_conductivity"],
self.vessel_density,
self.vessel_cp,
)
# Check if single-layer or composite (liner + shell) construction
if (
"liner_thermal_conductivity"
not in self.input["vessel"].keys()
):
# Single-layer wall construction
nn = 11 # number of nodes through wall thickness
z = np.linspace(0, self.thickness, nn)
self.z = z
# Create meshes for unwetted and wetted regions
mesh = tm.Mesh(z, tm.LinearElement)
mesh_w = tm.Mesh(z, tm.LinearElement)
# Material model with constant properties
cpeek = tm.isothermal_model(k, rho, cp)
cpeek_w = tm.isothermal_model(k, rho, cp)
# Initialize temperature profile on first time step
if type(T_profile) == type(int()) and T_profile == 0:
bc = [
{"T": self.T0},
{"T": self.Tamb},
]
domain = tm.Domain(mesh, [cpeek], bc)
domain.set_T(
(self.Tamb + self.T0) / 2 * np.ones(len(mesh.nodes))
)
solver = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded, T_profile = tm.solve_ht(domain, solver)
bc_w = [
{"T": self.T0},
{"T": self.Tamb},
]
domain_w = tm.Domain(mesh_w, [cpeek_w], bc_w)
domain_w.set_T(
(self.Tamb + self.T0)
/ 2
* np.ones(len(mesh_w.nodes))
)
solver_w = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded_w, T_profile_w = tm.solve_ht(
domain_w, solver_w
)
else:
# Boundary conditions: z=0 is outer wall, z=L is inner wall
bc = [
{
"q": self.q_outer[i]
# / (self.surf_area_outer - wetted_area_outer)
},
{
"q": -self.q_inner[i]
# / (self.surf_area_inner - wetted_area)
},
]
domain = tm.Domain(mesh, [cpeek], bc)
domain.set_T(T_profile[-1, :])
solver = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded, T_profile = tm.solve_ht(domain, solver)
# Wetted wall will be solved below if liquid is present
self.temp_profile.append(T_profile[-1, :])
self.T_outer_wall[i] = T_profile[-1, 0]
self.T_inner_wall[i] = T_profile[-1, -1]
# Only solve wetted wall if liquid is present
if self.liquid_level[i - 1] > 0 and wetted_area > 0:
bc_w = [
{
"q": self.q_outer_wetted[i]
}, # / wetted_area_outer},
{"q": -self.q_inner_wetted[i]}, # / wetted_area},
]
domain_w = tm.Domain(mesh_w, [cpeek_w], bc_w)
domain_w.set_T(T_profile_w[-1, :])
solver_w = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_w, T_profile_w = tm.solve_ht(domain_w, solver_w)
self.T_outer_wall_wetted[i] = T_profile_w[-1, 0]
self.T_inner_wall_wetted[i] = T_profile_w[-1, -1]
else:
# No liquid - wetted wall same as unwetted
self.T_outer_wall_wetted[i] = T_profile[-1, 0]
self.T_inner_wall_wetted[i] = T_profile[-1, -1]
else:
k_liner = self.input["vessel"]["liner_thermal_conductivity"]
rho_liner = self.input["vessel"]["liner_density"]
cp_liner = self.input["vessel"]["liner_heat_capacity"]
liner = tm.isothermal_model(k_liner, rho_liner, cp_liner)
shell = tm.isothermal_model(k, rho, cp)
liner_w = tm.isothermal_model(k_liner, rho_liner, cp_liner)
shell_w = tm.isothermal_model(k, rho, cp)
thk = self.input["vessel"]["thickness"] # thickness in m
nn = 11 # number of nodes
z_shell = np.linspace(0, thk, nn) # node locations
thk = self.input["vessel"]["liner_thickness"]
z_liner = np.linspace(-thk, 0, nn) # node locations
z2 = np.hstack((z_liner, z_shell[1:]))
self.z = z2
mesh2 = tm.Mesh(z2, tm.LinearElement)
mesh2_w = tm.Mesh(z2, tm.LinearElement)
for j, elem in enumerate(mesh2.elem):
if elem.nodes.mean() > 0.0:
mesh2.subdomain[j] = 1
mesh2_w.subdomain[j] = 1
if type(T_profile2) == type(int()) and T_profile2 == 0:
bc = [
{"T": self.T0},
{"T": self.Tamb},
]
domain2 = tm.Domain(mesh2, [liner, shell], bc)
domain2.set_T(
(self.Tamb + self.T0)
/ 2
* np.ones(len(mesh2.nodes))
)
solver2 = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded, T_profile2 = tm.solve_ht(domain2, solver2)
bc_w = [
{"T": self.T0},
{"T": self.Tamb},
]
domain2_w = tm.Domain(mesh2_w, [liner_w, shell_w], bc_w)
domain2_w.set_T(
(self.Tamb + self.T0)
/ 2
* np.ones(len(mesh2.nodes))
)
solver2_w = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded_w, T_profile2_w = tm.solve_ht(
domain2_w, solver2_w
)
else:
# Boundary conditions: z=-liner_thickness is inner, z=thickness is outer
bc = [
{
"q": -self.q_inner[i]
# / (self.surf_area_inner - wetted_area)
},
{
"q": self.q_outer[i]
# / (self.surf_area_outer - wetted_area_outer)
},
]
domain2 = tm.Domain(mesh2, [liner, shell], bc)
domain2.set_T(T_profile2[-1, :])
solver2 = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded, T_profile2 = tm.solve_ht(domain2, solver2)
bc_w = [
{"q": -self.q_inner_wetted[i]}, # / wetted_area},
{
"q": self.q_outer_wetted[i]
}, # / wetted_area_outer},
]
domain2_w = tm.Domain(mesh2_w, [liner_w, shell_w], bc_w)
domain2_w.set_T(T_profile2_w[-1, :])
solver2_w = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded_w, T_profile2_w = tm.solve_ht(
domain2_w, solver2_w
)
self.T_outer_wall[i] = T_profile2[-1, -1]
self.T_inner_wall[i] = T_profile2[-1, 0]
self.T_bonded_wall[i] = T_profile2[-1, (nn - 1)]
self.T_outer_wall_wetted[i] = T_profile2_w[-1, -1]
self.T_inner_wall_wetted[i] = T_profile2_w[-1, 0]
self.T_bonded_wall_wetted[i] = T_profile2_w[-1, (nn - 1)]
self.temp_profile.append(T_profile2[-1, :])
else:
# Lumped capacitance model (no 1D heat transfer)
self.T_inner_wall[i] = self.T_vessel[i]
self.T_outer_wall[i] = self.T_vessel[i]
self.T_inner_wall_wetted[i] = self.T_vessel_wetted[i]
self.T_outer_wall_wetted[i] = self.T_vessel_wetted[i]
# Calculate Biot number to validate lumped capacitance assumption
# Bi = h*L/k where L is characteristic length (wall thickness)
# Bi << 0.1: lumped model valid (uniform wall temperature)
# Bi > 0.1: thermal gradient significant, 1D model recommended
if "thickness" in self.input["vessel"]:
L_char = self.thickness
# Check for thermal_conductivity_biot first (for Biot calc only)
# Otherwise check thermal_conductivity (which triggers 1D model)
k_wall = None
if "thermal_conductivity_biot" in self.input["vessel"]:
k_wall = self.input["vessel"][
"thermal_conductivity_biot"
]
elif "thermal_conductivity" in self.input["vessel"]:
k_wall = self.input["vessel"]["thermal_conductivity"]
if k_wall is None:
self.Biot[i] = np.nan
self.Biot_wetted[i] = np.nan
else:
self.Biot[i] = hi * L_char / k_wall
self.Biot_wetted[i] = (
hiw * L_char / k_wall if hiw > 0 else 0.0
)
# Issue warning if Biot number exceeds threshold
if i == 1 or (i % 100 == 0 and self.Biot[i] > 0.1):
if self.Biot[i] > 0.1:
import warnings
warnings.warn(
f"t={self.time_array[i]:.1f}s: Biot number (unwetted) = {self.Biot[i]:.3f} > 0.1. "
f"Lumped capacitance model may be inaccurate. Consider using 1D heat transfer "
f"by specifying 'thermal_conductivity' in vessel properties.",
UserWarning,
)
if (
self.liquid_level[i - 1] > 0
and self.Biot_wetted[i] > 0.1
):
if i == 1 or i % 100 == 0:
import warnings
warnings.warn(
f"t={self.time_array[i]:.1f}s: Biot number (wetted) = {self.Biot_wetted[i]:.3f} > 0.1. "
f"Lumped capacitance model may be inaccurate.",
UserWarning,
)
else:
self.Biot[i] = np.nan
self.Biot_wetted[i] = np.nan
elif self.heat_method == "s-b":
if self.vessel_orientation == "horizontal":
L = self.diameter
else:
L = self.length
wetted_area = self.inner_vol.SA_from_h(self.liquid_level[i - 1])
if np.isnan(wetted_area):
wetted_area = 0
# Calculate outer wetted area for correct heat transfer on outer surface
# Add wall thickness to liquid level height since outer vessel reference
# is at bottom of wall, not bottom of inner surface
liquid_level_outer = self.liquid_level[i - 1] + self.thickness
wetted_area_outer = self.outer_vol.SA_from_h(liquid_level_outer)
if np.isnan(wetted_area_outer):
wetted_area_outer = 0
# Determine which wall temperature to use for internal heat transfer
# For 1D heat transfer: Use actual inner wall surface temperature
# For 0D model: Use bulk vessel temperature
if "thermal_conductivity" in self.input["vessel"].keys():
# 1D model: Use inner wall surface temperature
T_wall_inner = self.T_inner_wall[i - 1]
T_wall_inner_wetted = self.T_inner_wall_wetted[i - 1]
else:
# 0D model: Use bulk vessel temperature
T_wall_inner = self.T_vessel[i - 1]
T_wall_inner_wetted = self.T_vessel_wetted[i - 1]
# For NEM: use gas temperature for unwetted (gas-side) heat transfer
# For equilibrium: use bulk fluid temperature
if self.non_equilibrium and hasattr(self, 'T_gas'):
T_for_gas_side = self.T_gas[i - 1]
else:
T_for_gas_side = self.T_fluid[i - 1]
hi = tp.h_inner(
L,
T_for_gas_side,
T_wall_inner,
self.P[i - 1],
self.species,
)
self.h_inside[i] = hi
# NEM: Check liquid phase for boiling, use liquid properties
# Equilibrium: Use existing logic with equilibrium fluid
if self.non_equilibrium:
# For NEM, if liquid exists, use nucleate boiling correlation
# No need to check quality - liquid phase is always liquid
liquid_exists = self.m_liquid[i-1] > 1e-6
if liquid_exists:
# Use liquid temperature and liquid phase object
try:
self.transport_fluid_wet.update(
CP.PT_INPUTS, self.P[i - 1], self.T_liquid[i - 1]
)
except:
# If update fails, use saturation
self.transport_fluid_wet.update(CP.PQ_INPUTS, self.P[i - 1], 0.0)
# Update fluid_liquid to saturated state for h_inside_wetted
# (needed for surface tension and saturated properties)
self.fluid_liquid.update(CP.PQ_INPUTS, self.P[i - 1], 0.0)
hiw = tp.h_inside_wetted(
L,
T_wall_inner_wetted,
self.T_liquid[i - 1], # Use liquid temperature
self.transport_fluid_wet,
self.fluid_liquid, # Use liquid phase object at saturation
)
else:
# No liquid - use gas-side coefficient
hiw = hi
else:
# Equilibrium mode: use existing logic
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
self.transport_fluid_wet.update(
CP.PT_INPUTS, self.P[i - 1], self.T_fluid[i - 1]
)
hiw = tp.h_inside_wetted(
L,
T_wall_inner_wetted,
self.T_fluid[i - 1],
self.transport_fluid_wet,
self.fluid,
)
else:
hiw = hi
# For unwetted heat transfer, use gas temperature in NEM
if self.non_equilibrium and hasattr(self, 'T_gas'):
T_for_gas_side_htc = self.T_gas[i - 1]
else:
T_for_gas_side_htc = self.T_fluid[i - 1]
self.Q_inner[i] = self.scaling * (
(self.surf_area_inner - wetted_area)
* hi
* (T_wall_inner - T_for_gas_side_htc)
)
self.q_inner[i] = hi * (T_wall_inner - T_for_gas_side_htc)
# For wetted heat transfer, use liquid temperature in NEM
if self.non_equilibrium and hasattr(self, 'T_liquid'):
T_fluid_wet = self.T_liquid[i - 1]
else:
T_fluid_wet = self.T_fluid[i - 1]
self.Q_inner_wetted[i] = self.scaling * (
wetted_area * hiw * (T_wall_inner_wetted - T_fluid_wet)
)
self.q_inner_wetted[i] = hiw * (
T_wall_inner_wetted - T_fluid_wet
)
if np.isnan(self.Q_inner_wetted[i]):
self.Q_inner_wetted[i] = 0
# ================================================================
# FIRE HEAT LOAD CALCULATIONS (Stefan-Boltzmann Method)
# ================================================================
# Calculate external heat flux from fire using Stefan-Boltzmann
# equation accounting for radiative and convective heat transfer.
# Fire types: API 521 pool/jet fire, Scandpower pool/jet/peak fires
# Heat flux is temperature-dependent (higher vessel T → more re-radiation)
#
# For 1D heat transfer model: Use outer wall temperature (surface exposed to fire)
# For 0D model: Use bulk vessel temperature
# For two-phase systems: Calculate separate heat fluxes for
# wetted (liquid contact) and unwetted (gas contact) regions
# Determine which temperature to use for fire heat flux calculation
if "thermal_conductivity" in self.input["vessel"].keys():
# 1D model: Use actual outer wall surface temperature
T_fire_surface = self.T_outer_wall[i - 1]
T_fire_surface_wetted = self.T_outer_wall_wetted[i - 1]
else:
# 0D model: Use bulk vessel temperature
T_fire_surface = self.T_vessel[i - 1]
T_fire_surface_wetted = self.T_vessel_wetted[i - 1]
# Fire heats the outer surface - use outer surface area directly
self.Q_outer[i] = (
self.scaling
* fire.sb_fire(T_fire_surface, self.fire_type)
* (self.surf_area_outer - wetted_area_outer)
)
self.q_outer[i] = fire.sb_fire(T_fire_surface, self.fire_type)
self.Q_outer_wetted[i] = self.scaling * (
fire.sb_fire(T_fire_surface_wetted, self.fire_type)
* wetted_area_outer
)
self.q_outer_wetted[i] = fire.sb_fire(
T_fire_surface_wetted, self.fire_type
)
if np.isnan(self.Q_outer_wetted[i]):
self.Q_outer_wetted[i] = 0
# ================================================================
# 1-D TRANSIENT HEAT CONDUCTION (Fire Scenario with Detailed Thermal Model)
# ================================================================
# Solve transient heat conduction through vessel wall using
# finite element method (thermesh module) for fire scenarios.
# This section is activated when thermal_conductivity is specified.
# Otherwise, falls back to simple lumped capacitance model below.
if "thermal_conductivity" in self.input["vessel"].keys():
theta = 0.5 # Crank-Nicolson scheme (unconditionally stable, 2nd order)
dt = (
self.tstep / 10
) # Sub-step for thermal solver (finer time resolution)
k, rho, cp = (
self.input["vessel"]["thermal_conductivity"],
self.vessel_density,
self.vessel_cp,
)
# Check if single-layer or composite (liner + shell) construction
if (
"liner_thermal_conductivity"
not in self.input["vessel"].keys()
):
# Single-layer wall construction
nn = 11 # number of nodes through wall thickness
z = np.linspace(0, self.thickness, nn)
self.z = z
# Create meshes for unwetted and wetted regions
mesh = tm.Mesh(z, tm.LinearElement)
mesh_w = tm.Mesh(z, tm.LinearElement)
# Material model with constant properties
cpeek = tm.isothermal_model(k, rho, cp)
cpeek_w = tm.isothermal_model(k, rho, cp)
# Initialize temperature profile on first time step
if type(T_profile) == type(int()) and T_profile == 0:
bc = [
{"T": self.T0},
{"T": self.T0},
]
domain = tm.Domain(mesh, [cpeek], bc)
domain.set_T(self.T0 * np.ones(len(mesh.nodes)))
solver = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded, T_profile = tm.solve_ht(domain, solver)
bc_w = [
{"T": self.T0},
{"T": self.T0},
]
domain_w = tm.Domain(mesh_w, [cpeek_w], bc_w)
domain_w.set_T(self.T0 * np.ones(len(mesh_w.nodes)))
solver_w = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded_w, T_profile_w = tm.solve_ht(
domain_w, solver_w
)
else:
# Boundary conditions: z=0 is inner wall, z=L is outer wall
bc = [
{
"q": -self.q_inner[i]
# / (self.surf_area_inner - wetted_area)
},
{
"q": self.q_outer[i]
# / (self.surf_area_outer - wetted_area_outer)
},
]
domain = tm.Domain(mesh, [cpeek], bc)
domain.set_T(T_profile[-1, :])
solver = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded, T_profile = tm.solve_ht(domain, solver)
# Wetted wall boundary conditions will be set below
# in the conditional block that checks liquid_level
# Only solve wetted wall if liquid is present
# Check liquid_level to handle both gas-only and liquid-depleted cases
if self.liquid_level[i - 1] > 0 and wetted_area > 0:
bc_w = [
{"q": -self.q_inner_wetted[i]}, # / wetted_area},
{
"q": self.q_outer_wetted[i]
}, # , / wetted_area_outer},
]
domain_w = tm.Domain(mesh_w, [cpeek_w], bc_w)
domain_w.set_T(T_profile_w[-1, :])
solver_w = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_w, T_profile_w = tm.solve_ht(domain_w, solver_w)
self.T_inner_wall_wetted[i] = T_profile_w[-1, 0]
self.T_outer_wall_wetted[i] = T_profile_w[-1, -1]
self.T_vessel_wetted[i] = T_profile_w[-1, :].mean()
else:
# No liquid - wetted wall same as unwetted
self.T_inner_wall_wetted[i] = T_profile[-1, 0]
self.T_outer_wall_wetted[i] = T_profile[-1, -1]
self.T_vessel_wetted[i] = T_profile[-1, :].mean()
self.temp_profile.append(T_profile[-1, :])
self.T_inner_wall[i] = T_profile[-1, 0]
self.T_outer_wall[i] = T_profile[-1, -1]
# Update mean vessel temperature from wall temperatures
self.T_vessel[i] = T_profile[-1, :].mean()
else:
# Composite wall construction (liner + shell)
k_liner = self.input["vessel"]["liner_thermal_conductivity"]
rho_liner = self.input["vessel"]["liner_density"]
cp_liner = self.input["vessel"]["liner_heat_capacity"]
liner = tm.isothermal_model(k_liner, rho_liner, cp_liner)
shell = tm.isothermal_model(k, rho, cp)
liner_w = tm.isothermal_model(k_liner, rho_liner, cp_liner)
shell_w = tm.isothermal_model(k, rho, cp)
thk = self.input["vessel"]["thickness"] # thickness in m
nn = 11 # number of nodes
z_shell = np.linspace(0, thk, nn) # node locations
thk = self.input["vessel"]["liner_thickness"]
z_liner = np.linspace(-thk, 0, nn) # node locations
z2 = np.hstack((z_liner, z_shell[1:]))
self.z = z2
mesh2 = tm.Mesh(z2, tm.LinearElement)
mesh2_w = tm.Mesh(z2, tm.LinearElement)
for j, elem in enumerate(mesh2.elem):
if elem.nodes.mean() > 0.0:
mesh2.subdomain[j] = 1
mesh2_w.subdomain[j] = 1
if type(T_profile2) == type(int()) and T_profile2 == 0:
bc = [
{"T": self.T0},
{"T": self.T0},
]
domain2 = tm.Domain(mesh2, [liner, shell], bc)
domain2.set_T(self.T0 * np.ones(len(mesh2.nodes)))
solver2 = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded, T_profile2 = tm.solve_ht(domain2, solver2)
bc_w = [
{"T": self.T0},
{"T": self.T0},
]
domain2_w = tm.Domain(mesh2_w, [liner_w, shell_w], bc_w)
domain2_w.set_T(self.T0 * np.ones(len(mesh2.nodes)))
solver2_w = {
"dt": 100,
"t_end": 10000,
"theta": theta,
}
t_bonded_w, T_profile2_w = tm.solve_ht(
domain2_w, solver2_w
)
else:
# Boundary conditions: z=-liner_thickness is inner, z=thickness is outer
bc = [
{
"q": -self.q_inner[i]
# / (self.surf_area_inner - wetted_area)
},
{
"q": self.q_outer[i]
# / (self.surf_area_outer - wetted_area_outer)
},
]
domain2 = tm.Domain(mesh2, [liner, shell], bc)
domain2.set_T(T_profile2[-1, :])
solver2 = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded, T_profile2 = tm.solve_ht(domain2, solver2)
# Only solve wetted wall if liquid is present
# Check liquid_level to handle both gas-only and liquid-depleted cases
if self.liquid_level[i - 1] > 0 and wetted_area > 0:
bc_w = [
{
"q": -self.q_inner_wetted[i]
}, # / wetted_area},
{
"q": self.q_outer_wetted[i]
# / wetted_area_outer
},
]
domain2_w = tm.Domain(
mesh2_w, [liner_w, shell_w], bc_w
)
domain2_w.set_T(T_profile2_w[-1, :])
solver2_w = {
"dt": dt,
"t_end": self.tstep,
"theta": theta,
}
t_bonded_w, T_profile2_w = tm.solve_ht(
domain2_w, solver2_w
)
self.T_inner_wall_wetted[i] = T_profile2_w[-1, 0]
self.T_outer_wall_wetted[i] = T_profile2_w[-1, -1]
self.T_bonded_wall_wetted[i] = T_profile2_w[
-1, (nn - 1)
]
self.T_vessel_wetted[i] = T_profile2_w[-1, :].mean()
else:
# No liquid - wetted wall same as unwetted
self.T_inner_wall_wetted[i] = T_profile2[-1, 0]
self.T_outer_wall_wetted[i] = T_profile2[-1, -1]
self.T_bonded_wall_wetted[i] = T_profile2[
-1, (nn - 1)
]
self.T_vessel_wetted[i] = T_profile2[-1, :].mean()
self.T_inner_wall[i] = T_profile2[-1, 0]
self.T_outer_wall[i] = T_profile2[-1, -1]
self.T_bonded_wall[i] = T_profile2[-1, (nn - 1)]
self.temp_profile.append(T_profile2[-1, :])
# Update mean vessel temperature from wall temperatures
self.T_vessel[i] = T_profile2[-1, :].mean()
else:
# ================================================================
# SIMPLE LUMPED CAPACITANCE MODEL (No thermal gradient)
# ================================================================
# When thermal_conductivity is not specified, use simple 0D model
# with uniform vessel wall temperature (no spatial gradient)
self.T_vessel[i] = self.T_vessel[i - 1] + (
(self.Q_outer[i] - self.Q_inner[i]) / self.scaling
) * self.tstep / (
self.vessel_cp
* self.vessel_density
* self.vol_solid
* (self.inner_vol.A - wetted_area)
/ self.inner_vol.A
)
if self.liquid_level[i - 1] > 0:
self.T_vessel_wetted[i] = self.T_vessel_wetted[i - 1] + (
(self.Q_outer_wetted[i] - self.Q_inner_wetted[i])
/ self.scaling
) * self.tstep / (
self.vessel_cp
* self.vessel_density
* self.vol_solid
* wetted_area
/ self.inner_vol.A
)
else:
# Hack to heat up previous liquid wetted surface
# Fire heats outer surface - use outer surface area directly
self.T_vessel_wetted[i] = self.T_vessel_wetted[i - 1] + (
fire.sb_fire(
self.T_vessel_wetted[i - 1], self.fire_type
)
* (self.surf_area_outer - wetted_area_outer)
- (self.surf_area_inner - wetted_area)
* hi
* (self.T_vessel_wetted[i - 1] - self.T_fluid[i - 1])
) * self.tstep / (
self.vessel_cp
* self.vessel_density
* self.vol_solid
* (self.inner_vol.A - wetted_area)
/ self.inner_vol.A
)
if np.isnan(self.T_vessel_wetted[i]):
self.T_vessel_wetted[i] = self.T_vessel[i]
# Calculate Biot number to validate lumped capacitance assumption
# Bi = h*L/k where L is characteristic length (wall thickness)
# Bi << 0.1: lumped model valid (uniform wall temperature)
# Bi > 0.1: thermal gradient significant, 1D model recommended
if "thickness" in self.input["vessel"]:
L_char = (
self.thickness
) # Characteristic length = wall thickness
# Use inner heat transfer coefficient (typically controls)
# and vessel thermal conductivity if available
# Check for thermal_conductivity_biot first (for Biot calc only)
# Otherwise check thermal_conductivity (which triggers 1D model)
k_wall = None
if "thermal_conductivity_biot" in self.input["vessel"]:
k_wall = self.input["vessel"][
"thermal_conductivity_biot"
]
elif "thermal_conductivity" in self.input["vessel"]:
k_wall = self.input["vessel"]["thermal_conductivity"]
if k_wall is None:
# Lumped model without k specified - can't calc Bi, set to NaN
self.Biot[i] = np.nan
self.Biot_wetted[i] = np.nan
else:
self.Biot[i] = hi * L_char / k_wall
self.Biot_wetted[i] = (
hiw * L_char / k_wall if hiw > 0 else 0.0
)
# Issue warning if Biot number exceeds threshold
if i == 1 or (i % 100 == 0 and self.Biot[i] > 0.1):
if self.Biot[i] > 0.1:
import warnings
warnings.warn(
f"t={self.time_array[i]:.1f}s: Biot number (unwetted) = {self.Biot[i]:.3f} > 0.1. "
f"Lumped capacitance model may be inaccurate. Consider using 1D heat transfer "
f"by specifying 'thermal_conductivity' in vessel properties.",
UserWarning,
)
if (
self.liquid_level[i - 1] > 0
and self.Biot_wetted[i] > 0.1
):
if i == 1 or i % 100 == 0:
import warnings
warnings.warn(
f"t={self.time_array[i]:.1f}s: Biot number (wetted) = {self.Biot_wetted[i]:.3f} > 0.1. "
f"Lumped capacitance model may be inaccurate.",
UserWarning,
)
else:
self.Biot[i] = np.nan
self.Biot_wetted[i] = np.nan
self.T_inner_wall[i] = self.T_vessel[i]
self.T_outer_wall[i] = self.T_vessel[i]
self.T_inner_wall_wetted[i] = self.T_vessel_wetted[i]
self.T_outer_wall_wetted[i] = self.T_vessel_wetted[i]
elif self.heat_method == "specified_U":
self.Q_inner[i] = (
self.surf_area_outer
* self.Ufix
* (self.Tamb - self.T_fluid[i - 1])
)
self.T_vessel[i] = self.T_vessel[0]
elif self.heat_method == "specified_Q":
self.Q_inner[i] = self.Qfix
self.T_vessel[i] = self.T_vessel[0]
else:
self.Q_inner[i] = 0.0
self.T_vessel[i] = self.T_vessel[0]
# NMOL = self.mass_fluid[i - 1] / self.MW
# NMOL_ADD = (self.mass_fluid[i] - self.mass_fluid[i - 1]) / self.MW
# New
U_start = self.U_mass[i - 1] * self.mass_fluid[i - 1]
# Smooting finction for very early times /numerical trick
# Might not be necessary.
x = 1 - math.exp(-1 * self.time_array[i]) ** 0.66
# Finding the inlet/outlet enthalpy rate for the energy balance
if input["valve"]["flow"] == "filling":
# h_in = self.fluid.hmass()
h_in = x * self.res_fluid.hmass() + (1 - x) * self.res_fluid.umass()
else:
# h_in = self.fluid.hmass()
if self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
h_in = self.fluid.saturated_vapor_keyed_output(CP.iHmass)
else:
h_in = self.fluid.hmass()
if i > 1:
P1 = self.P[i - 2]
else:
P1 = self.P[i - 1]
U_end = (
U_start
- self.tstep * self.mass_rate[i - 1] * h_in
+ self.tstep * self.Q_inner[i]
+ self.tstep * self.Q_inner_wetted[i]
)
self.U_mass[i] = U_end / self.mass_fluid[i]
# ====================================================================
# NON-EQUILIBRIUM MODEL (NEM) ENERGY BALANCE
# ====================================================================
# For non-equilibrium model, solve separate energy balances for gas and liquid phases
if self.non_equilibrium:
# First, update masses with valve flow only (no phase transfer yet)
# Assume discharge/filling happens through gas phase
if input["valve"]["flow"] == "discharge":
# Mass leaves through gas phase
self.m_gas[i] = self.m_gas[i-1] - self.mass_rate[i-1] * self.tstep
self.m_liquid[i] = self.m_liquid[i-1]
else: # filling
# Mass enters through gas phase
self.m_gas[i] = self.m_gas[i-1] - self.mass_rate[i-1] * self.tstep
self.m_liquid[i] = self.m_liquid[i-1]
# Prevent negative masses
if self.m_gas[i] < 0:
self.m_gas[i] = 0.0
if self.m_liquid[i] < 0:
self.m_liquid[i] = 0.0
# ====================================================================
# PHASE TRANSFER (Part of Mass Balance)
# ====================================================================
# Check if previous timestep had two-phase conditions in either phase
# Transfer mass and energy between phases based on vapor quality
# This happens BEFORE energy balance, as part of mass redistribution
dm_evap_mass = 0.0 # Mass evaporated (liquid → gas)
dm_cond_mass = 0.0 # Mass condensed (gas → liquid)
E_evap = 0.0 # Energy transferred with evaporation
E_cond = 0.0 # Energy transferred with condensation
# Check liquid phase from previous timestep
if self.m_liquid[i-1] > 1e-6:
try:
# Check if liquid was in two-phase region at previous state
self.fluid_liquid.update(CP.DmassUmass_INPUTS, self.rho_liquid[i-1], self.U_liquid[i-1])
quality_liquid = self.fluid_liquid.Q()
if quality_liquid > 0 and quality_liquid < 1:
# Liquid is in two-phase region - some should evaporate
relax_factor = 0.8 # Transfer 80% per timestep
dm_evap_mass = quality_liquid * self.m_liquid[i-1] * relax_factor
# Safety limit: evaporation rate limited by available heat
# Maximum dm/dt from heat input: Q / h_fg
h_liq_sat_check = self.fluid_liquid.saturated_liquid_keyed_output(CP.iHmass)
h_vap_sat_check = self.fluid_liquid.saturated_vapor_keyed_output(CP.iHmass)
h_fg = h_vap_sat_check - h_liq_sat_check
# Estimate available heat for evaporation (use PREVIOUS timestep)
Q_available = max(self.Q_inner_wetted[i-1], 0.0) * self.tstep # J
dm_evap_max = Q_available / h_fg if h_fg > 0 else dm_evap_mass
# Limit to 50x heat-based rate to allow for stored energy
dm_evap_mass = min(dm_evap_mass, dm_evap_max * 50.0)
# Transfer mass
self.m_liquid[i] -= dm_evap_mass
self.m_gas[i] += dm_evap_mass
# Energy transferred with phase change
# Evaporated mass leaves liquid as VAPOR (phase change occurs)
# The evaporated liquid becomes saturated vapor
# Use (h+u)/2 compromise between enthalpy and internal energy
h_vap_sat = self.fluid_liquid.saturated_vapor_keyed_output(CP.iHmass)
u_vap_sat = self.fluid_liquid.saturated_vapor_keyed_output(CP.iUmass)
e_vap_sat = (h_vap_sat + u_vap_sat) / 2.0
# E_evap = energy carried by evaporated mass
# This energy is SUBTRACTED from liquid and ADDED to gas
E_evap = dm_evap_mass * e_vap_sat
except:
pass # Not in two-phase, no transfer
# Check gas phase from previous timestep (only if no evaporation)
if self.m_gas[i-1] > 1e-6 and dm_evap_mass == 0:
try:
# Check if gas was in two-phase region at previous state
self.fluid_gas.update(CP.DmassUmass_INPUTS, self.rho_gas[i-1], self.U_gas[i-1])
quality_gas = self.fluid_gas.Q()
if quality_gas > 0 and quality_gas < 1:
# Gas is in two-phase region - some should condense
relax_factor = 0.8 # Transfer 80% per timestep
dm_cond_mass = (1.0 - quality_gas) * self.m_gas[i-1] * relax_factor
# Safety limit: condensation rate limited by heat removal
h_liq_sat_check = self.fluid_gas.saturated_liquid_keyed_output(CP.iHmass)
h_vap_sat_check = self.fluid_gas.saturated_vapor_keyed_output(CP.iHmass)
h_fg = h_vap_sat_check - h_liq_sat_check
# Estimate heat removal rate (use PREVIOUS timestep)
Q_removal = max(-self.Q_inner[i-1], 0.0) * self.tstep # J (positive value)
dm_cond_max = Q_removal / h_fg if h_fg > 0 else dm_cond_mass
# Limit to 50x heat-based rate
dm_cond_mass = min(dm_cond_mass, dm_cond_max * 50.0)
# Transfer mass
self.m_gas[i] -= dm_cond_mass
self.m_liquid[i] += dm_cond_mass
# Energy transferred with phase change
# Condensed mass leaves gas as LIQUID (phase change occurs)
# The condensed vapor becomes saturated liquid
# Use (h+u)/2 compromise between enthalpy and internal energy
h_liq_sat = self.fluid_gas.saturated_liquid_keyed_output(CP.iHmass)
u_liq_sat = self.fluid_gas.saturated_liquid_keyed_output(CP.iUmass)
e_liq_sat = (h_liq_sat + u_liq_sat) / 2.0
# E_cond = energy carried by condensed mass
# This energy is SUBTRACTED from gas and ADDED to liquid
E_cond = dm_cond_mass * e_liq_sat
except:
pass # Not in two-phase, no transfer
# Store phase transfer rate for output
net_transfer = dm_cond_mass - dm_evap_mass # Positive = condensation
self.mdot_phase_transfer[i] = net_transfer / self.tstep
# Calculate gas-liquid interfacial heat transfer
# Q_gl = h_gl * A_interface * (T_gas - T_liquid)
# Positive = heat from gas to liquid
if self.m_liquid[i-1] > 1e-6 and self.m_gas[i-1] > 1e-6:
# Calculate interface area from liquid level using TANK geometry
V_liquid_prev = self.m_liquid[i-1] / self.rho_liquid[i-1]
ll_prev = self.inner_vol.h_from_V(V_liquid_prev)
# Use fluids.TANK.A_cross_sectional for exact interface area
# Handles horizontal/vertical, all head types (hemispherical, F&D, etc.)
A_interface = self.inner_vol.A_cross_sectional(ll_prev)
# Heat transfer coefficient (W/m²K)
# Can be specified in input file, calculated, or use default
if "h_gas_liquid" in input["calculation"]:
h_gl_input = input["calculation"]["h_gas_liquid"]
if isinstance(h_gl_input, str) and h_gl_input.lower() in ("calc", "calc_two_sided"):
# Calculate h_gl using natural convection correlation
L_char = self.inner_vol.D
try:
if h_gl_input.lower() == "calc_two_sided":
h_gl = tp.h_gas_liquid_interface_two_sided(
self.T_gas[i-1],
self.T_liquid[i-1],
self.P[i-1],
L_char,
self.fluid_gas,
self.fluid_liquid
)
else:
h_gl = tp.h_gas_liquid_interface(
self.T_gas[i-1],
self.T_liquid[i-1],
self.P[i-1],
L_char,
self.fluid_gas,
self.fluid_liquid
)
except Exception as e:
if i == 1 or (i < 100 and i % 50 == 0):
import warnings
warnings.warn(f"h_gl correlation failed at t={self.time_array[i]:.1f}s: {str(e)[:80]}, using fallback h_gl=100 W/m²K")
h_gl = 100.0
else:
# User-specified fixed value
h_gl = float(h_gl_input)
else:
h_gl = 50.0 # Default value
# Store h_gl in array for diagnostics
self.h_gas_liquid[i] = h_gl
# Heat transfer rate (W)
Q_gas_liquid = h_gl * A_interface * (self.T_gas[i-1] - self.T_liquid[i-1])
else:
Q_gas_liquid = 0.0
self.h_gas_liquid[i] = 0.0
# ====================================================================
# ENERGY BALANCE
# ====================================================================
# Apply heat transfer and flow work to each phase
# Phase transfer energy already accounted for in mass balance
# Gas phase energy balance:
# dU_gas = Q_wall_gas - Q_gas_liquid + h_valve*dm_valve + E_phase_transfer
U_gas_start = self.U_gas[i-1] * self.m_gas[i-1]
# Enthalpy for valve flow
if input["valve"]["flow"] == "filling":
h_valve = self.res_fluid.hmass()
U_gas_tentative = (U_gas_start
- self.tstep * self.mass_rate[i-1] * h_valve
+ self.tstep * self.Q_inner[i]
- self.tstep * Q_gas_liquid
+ E_evap - E_cond) # Phase transfer energy
else: # discharge
if self.m_gas[i-1] > 1e-6:
# Use DmassUmass to get enthalpy (avoids PT issues at saturation)
# Calculate discharge enthalpy - fail if thermodynamic state is invalid
self.fluid_gas.update(CP.DmassUmass_INPUTS, self.rho_gas[i-1], self.U_gas[i-1])
h_valve = self.fluid_gas.hmass()
else:
h_valve = 0.0
U_gas_tentative = (U_gas_start
- self.tstep * self.mass_rate[i-1] * h_valve
+ self.tstep * self.Q_inner[i]
- self.tstep * Q_gas_liquid
+ E_evap - E_cond) # Phase transfer energy
# Liquid phase energy balance:
# dU_liquid = Q_wall_liquid + Q_gas_liquid - E_phase_transfer
U_liquid_start = self.U_liquid[i-1] * self.m_liquid[i-1]
U_liquid_tentative = (U_liquid_start
+ self.tstep * self.Q_inner_wetted[i]
+ self.tstep * Q_gas_liquid
- E_evap + E_cond) # Phase transfer energy (opposite sign)
# Final internal energies for P-solver
U_gas_end = U_gas_tentative
U_liquid_end = U_liquid_tentative
# Safety check: if gas mass is very small, switch to single-phase liquid
if self.m_gas[i] < 1e-3: # Less than 1 gram of gas
# Essentially single-phase liquid - add remaining gas energy to liquid
U_liquid_end += U_gas_end
U_gas_end = 0.0
self.m_gas[i] = 0.0
# Safety check: if liquid mass is very small, switch to single-phase gas
if self.m_liquid[i] < 1e-3: # Less than 1 gram of liquid
# Essentially single-phase gas - add remaining liquid energy to gas
U_gas_end += U_liquid_end
U_liquid_end = 0.0
self.m_liquid[i] = 0.0
# Solve for pressure P such that volume constraint is satisfied
# Given P and U, CoolProp directly gives T and ρ (no nested iteration needed!)
from scipy.optimize import brentq, brenth, ridder, bisect, minimize, newton
def volume_residual(P):
"""
For given pressure P, update phases with (P, U) and check volume constraint.
Returns (V_gas + V_liquid - V_vessel) / V_vessel
"""
try:
U_gas_specific = U_gas_end / self.m_gas[i]
U_liquid_specific = U_liquid_end / self.m_liquid[i]
# Update gas with (P, U) → CoolProp gives ρ and T directly
self.fluid_gas.update(CP.PUmass_INPUTS, P, U_gas_specific)
rho_gas = self.fluid_gas.rhomass()
# Update liquid with (P, U) → CoolProp gives ρ and T directly
self.fluid_liquid.update(CP.PUmass_INPUTS, P, U_liquid_specific)
rho_liquid = self.fluid_liquid.rhomass()
# Calculate volumes
V_gas = self.m_gas[i] / rho_gas
V_liquid = self.m_liquid[i] / rho_liquid
V_total = V_gas + V_liquid
# Return normalized volume residual
return ((V_total - self.vol) / self.vol)
except Exception:
# CoolProp PUmass_INPUTS can fail near the critical point.
# Return a sign-consistent residual based on physical behavior:
# low P → expanded gas → V_total > V_vessel → positive
# high P → compressed gas → V_total < V_vessel → negative
# Using +1.0 for all failures creates false sign changes
# that mislead bracket-based solvers (brentq/bisect).
if P > self.P[i-1]:
return -10.0 # High P side: compressed → negative
else:
return 10.0 # Low P side: expanded → positive
try:
# Check if we're essentially single-phase
if self.m_gas[i] < 1e-6 and self.m_liquid[i] > 1e-6:
# Single-phase liquid - use DU update
U_total_specific = U_liquid_end / self.m_liquid[i]
rho_total = self.mass_fluid[i] / self.vol
self.fluid.update(CP.DmassUmass_INPUTS, rho_total, U_total_specific)
P_solution = self.fluid.p()
elif self.m_liquid[i] < 1e-6 and self.m_gas[i] > 1e-6:
# Single-phase gas - use DU update
U_total_specific = U_gas_end / self.m_gas[i]
rho_total = self.mass_fluid[i] / self.vol
self.fluid.update(CP.DmassUmass_INPUTS, rho_total, U_total_specific)
P_solution = self.fluid.p()
else:
# Two-phase - solve for pressure equilibrium
# Two-phase VLE physically cannot exist above P_crit,
# so capping bounds at P_crit*0.95 is correct here.
# (When liquid depletes, code switches to single-phase
# DU path above, which has no pressure cap.)
P_crit = self.fluid_gas.p_critical()
P_min = max(self.P[i-1] * 0.5, 1e5) # At least 1 bar
P_max = min(self.P[i-1] * 2.0, P_crit * 0.95)
# Solve for pressure that gives correct volume
# Use fallback solver chain for robustness
P_solution = None
solver_method = None
try:
# Method 1: Try brentq (fastest)
P_solution = brentq(volume_residual, P_min, P_max, xtol=1e-5, maxiter=100)
solver_method = "brentq"
except ValueError:
# Root not bracketed - try ridder (more robust)
try:
P_solution = ridder(volume_residual, P_min, P_max, xtol=1e-5, maxiter=100)
solver_method = "ridder"
except ValueError:
# Still not bracketed - expand bounds and try bisect
P_min_expanded = max(P_min * 0.2, 1e5) # Expand to 0.2x
P_max_expanded = min(P_max * 5.0, P_crit * 0.95)
try:
P_solution = bisect(volume_residual, P_min_expanded, P_max_expanded, xtol=1e-5, maxiter=200)
solver_method = "bisect_expanded"
except ValueError:
# Last resort: check if previous pressure gives acceptable residual
res_prev = volume_residual(self.P[i-1])
if abs(res_prev) < 0.05: # Within 5% of target volume
P_solution = self.P[i-1]
solver_method = "prev_pressure"
else:
# Give up - cannot find solution
raise ValueError(
f"All solvers failed. Residual at P_prev={self.P[i-1]/1e5:.2f} bar: {res_prev:.4f}"
)
# Store solved pressure
self.P[i] = P_solution
# Update gas state with solved pressure
if self.m_gas[i] > 1e-6:
U_gas_specific = U_gas_end / self.m_gas[i]
self.fluid_gas.update(CP.PUmass_INPUTS, P_solution, U_gas_specific)
self.rho_gas[i] = self.fluid_gas.rhomass()
self.U_gas[i] = self.fluid_gas.umass()
self.T_gas[i] = self.fluid_gas.T()
else:
self.rho_gas[i] = 0.0
self.U_gas[i] = 0.0
self.T_gas[i] = self.T_gas[i-1]
# Update liquid state with solved pressure
if self.m_liquid[i] > 1e-6:
U_liquid_specific = U_liquid_end / self.m_liquid[i]
self.fluid_liquid.update(CP.PUmass_INPUTS, P_solution, U_liquid_specific)
self.rho_liquid[i] = self.fluid_liquid.rhomass()
self.U_liquid[i] = self.fluid_liquid.umass()
self.T_liquid[i] = self.fluid_liquid.T()
else:
self.rho_liquid[i] = 0.0
self.U_liquid[i] = 0.0
self.T_liquid[i] = self.T_liquid[i-1]
# Update combined state for compatibility
self.T_fluid[i] = self.T_gas[i] if self.m_gas[i] > 1e-6 else self.T_liquid[i]
self.rho[i] = self.mass_fluid[i] / self.vol
# Update main fluid object for compatibility
total_U = (U_gas_end + U_liquid_end) / self.mass_fluid[i]
self.fluid.update(CP.DmassUmass_INPUTS, self.rho[i], total_U)
# Calculate liquid level for NEM
if self.m_liquid[i] > 1e-6 and self.rho_liquid[i] > 1e-6:
V_liquid = self.m_liquid[i] / self.rho_liquid[i]
self.liquid_level[i] = self.inner_vol.h_from_V(V_liquid)
else:
self.liquid_level[i] = 0.0
# Phase transfer already handled in mass balance step
# No post-solver adjustments needed
except Exception as e:
# Fail hard - do not allow pressure equilibrium violations
raise RuntimeError(
f"NEM solver failed at t={self.time_array[i]:.2f}s: {e}\n"
f" P_guess range: [{P_min:.2f}, {P_max:.2f}] Pa\n"
f" m_gas: {self.m_gas[i]:.3f} kg, m_liquid: {self.m_liquid[i]:.3f} kg\n"
f" U_gas_end: {U_gas_end:.1f} J, U_liquid_end: {U_liquid_end:.1f} J\n"
f"Solver must maintain strict pressure equilibrium. Check:\n"
f" 1. Timestep may be too large (current: {self.tstep}s)\n"
f" 2. Energy balance may have unphysical values\n"
f" 3. Phase transfer rate may be too aggressive"
) from e
# Not pretty if-statement and a hack for fire relief area estimation. Most cases go directly to the first ...else... clause
elif input["valve"]["type"] == "relief":
if self.Pset <= self.P[i - 1]:
if self.liquid_level[i - 1] > 0:
self.fluid.update(
CP.HmassP_INPUTS,
self.fluid.hmass()
+ self.tstep
* (self.Q_inner[i] + self.Q_inner_wetted[i])
/ self.mass_fluid[i],
self.Pset,
)
total_mass = self.fluid.rhomass() * self.vol
Hvap = self.fluid.saturated_vapor_keyed_output(
CP.iHmass
) - self.fluid.saturated_liquid_keyed_output(CP.iHmass)
BOG = (
(self.Q_inner_wetted[i] + self.Q_inner[i])
/ (Hvap)
* 3600
* 24
/ self.mass_fluid[i]
* 100
)
# print("BOG rate at relief (%%/day): ", BOG)
self.mass_rate[i] = (
self.Q_inner_wetted[i] + self.Q_inner[i]
) / Hvap
P1 = self.Pset
self.P[i] = P1
self.T_fluid[i] = self.fluid.T()
Z = self.fluid.saturated_vapor_keyed_output(CP.iZ)
cp0molar = self.fluid.saturated_vapor_keyed_output(
CP.iCp0molar
)
self.relief_area[i] = fluids.API520_A_g(
self.mass_rate[i],
self.fluid.T(),
Z,
self.MW * 1000,
cp0molar / (cp0molar - 8.314),
P1,
self.p_back,
0.975,
1,
1,
)
else:
T1 = self.PHproblem(
h_in
+ self.tstep * self.Q_inner[i] / self.mass_fluid[i],
self.Pset,
Tguess=self.T_fluid[i - 1] + 5,
relief=True,
)
self.T_fluid[i] = T1
P1 = self.Pset
self.P[i] = P1
self.T_fluid[i] = T1
self.fluid.update(CP.PT_INPUTS, self.P[i], T1)
self.mass_rate[i] = (
(
(1 / self.fluid.rhomass() - 1 / self.rho[i])
* self.mass_fluid[i]
)
* self.rho[i]
/ self.tstep
)
self.relief_area[i] = fluids.API520_A_g(
self.mass_rate[i],
T1,
self.fluid.compressibility_factor(),
self.MW * 1000,
self.fluid.cp0molar() / (self.fluid.cp0molar() - 8.314),
P1,
self.p_back,
0.975,
1,
1,
)
else:
if self.liquid_level[i - 1] > 0:
self.mass_rate[i] = 0
# self.fluid.update(CP.PQ_INPUTS, self.P[i], Q)
self.fluid.update(
CP.DmassUmass_INPUTS,
self.rho[i],
U_end / self.mass_fluid[i],
)
self.T_fluid[i] = self.fluid.T()
self.P[i] = self.fluid.p()
else:
self.mass_rate[i] = 0
P1, T1, self.U_res[i] = self.UDproblem(
U_end / self.mass_fluid[i],
self.rho[i],
self.P[i - 1],
self.T_fluid[i - 1],
)
self.P[i] = P1
self.T_fluid[i] = T1
self.fluid.update(CP.PT_INPUTS, self.P[i], self.T_fluid[i])
else:
P1, T1, self.U_res[i] = self.UDproblem(
U_end / self.mass_fluid[i],
self.rho[i],
self.P[i - 1],
self.T_fluid[i - 1],
)
self.P[i] = P1
self.T_fluid[i] = T1
if len(self.molefracs) == 1 and self.molefracs[0] == 1.0:
self.fluid.update(
CP.DmassUmass_INPUTS,
self.rho[i],
self.U_mass[i],
)
else:
try:
self.fluid.update(CP.PT_INPUTS, self.P[i], self.T_fluid[i])
except:
if self.fluid.Q() < 0:
self.fluid.update(CP.PQ_INPUTS, self.P[i], 1)
else:
self.fluid.update(
CP.PQ_INPUTS, self.P[i], self.fluid.Q()
)
if (
self.input["valve"]["flow"] == "discharge"
and self.fluid.Q() < 1
and self.fluid.Q() >= 0
):
self.res_fluid.update(CP.PQ_INPUTS, self.P[i], 1.0)
else:
raise NameError("Unknown calculation method: " + self.method)
Q = self.fluid.Q()
if Q >= 0 and Q <= 1:
self.vapour_mole_fraction[i] = Q
else:
self.vapour_mole_fraction[i] = 1
self.H_mass[i] = self.fluid.hmass()
self.S_mass[i] = self.fluid.smass()
self.U_mass[i] = self.fluid.umass()
self.liquid_level[i] = self.calc_liquid_level()
# Calculating vent temperature (adiabatic) only for discharge problem
if self.input["valve"]["flow"] == "discharge":
if "&" in self.species:
self.T_vent[i] = self.PHproblem(
self.H_mass[i], self.p_back, self.vent_fluid.T()
)
else:
try:
self.T_vent[i] = PropsSI(
"T", "H", self.H_mass[i], "P", self.p_back, self.species
)
except:
self.T_vent[i] = self.vent_fluid.T()
# For NEM: use gas phase properties for discharge calculations
# For equilibrium: use main fluid object
if self.non_equilibrium and self.m_gas[i] > 1e-6:
# NEM: check if gas phase is two-phase or superheated
Q_gas = self.fluid_gas.Q()
if Q_gas >= 0 and Q_gas <= 1:
# Two-phase: use saturated vapor properties
cpcv = self.fluid_gas.saturated_vapor_keyed_output(CP.iCpmolar) / (
self.fluid_gas.saturated_vapor_keyed_output(CP.iCpmolar) - 8.314
)
Z = self.fluid_gas.saturated_vapor_keyed_output(CP.iZ)
else:
# Superheated gas: use actual gas phase properties
cpcv = self.fluid_gas.cp0molar() / (self.fluid_gas.cp0molar() - 8.314)
Z = self.fluid_gas.compressibility_factor()
elif self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
# Equilibrium two-phase: use saturated vapor properties
cpcv = self.fluid.saturated_vapor_keyed_output(CP.iCpmolar) / (
self.fluid.saturated_vapor_keyed_output(CP.iCpmolar) - 8.314
)
Z = self.fluid.saturated_vapor_keyed_output(CP.iZ)
else:
# Single phase gas: use actual fluid properties
cpcv = self.fluid.cp0molar() / (self.fluid.cp0molar() - 8.314)
Z = self.fluid.compressibility_factor()
# ====================================================================
# MASS FLOW RATE CALCULATIONS (End of Time Step)
# ====================================================================
# Calculate mass flow rate for next time step based on valve type
# and current thermodynamic state. Uses current pressure, temperature,
# density to determine compressible flow through valve/orifice.
#
# Valve types:
# - orifice: Yellow Book compressible flow equation (critical/subcritical)
# - controlvalve: Control valve Cv characteristic with time-dependent opening
# - psv: Relief valve with API 520/521 sizing equations
# - mdot: Constant or time-varying mass flow rate (already set)
#
# For two-phase systems: Use saturated vapor properties for mass flow calc
# Note: "relief" valve type handled separately (not here)
if input["valve"]["type"] == "orifice":
if input["valve"]["flow"] == "filling":
k = self.res_fluid.cp0molar() / (self.res_fluid.cp0molar() - 8.314)
self.mass_rate[i] = -tp.gas_release_rate(
self.p_back,
self.P[i],
self.res_fluid.rhomass(),
k,
self.CD,
self.D_orifice**2 / 4 * math.pi,
)
else:
# For NEM: use gas phase density (actual state, not always saturated)
# For equilibrium two-phase: use saturated vapor density
if self.non_equilibrium and self.m_gas[i] > 1e-6:
Q_gas = self.fluid_gas.Q()
if Q_gas >= 0 and Q_gas <= 1:
rho = self.fluid_gas.saturated_vapor_keyed_output(CP.iDmass)
else:
rho = self.fluid_gas.rhomass()
elif self.fluid.Q() >= 0 and self.fluid.Q() <= 1:
rho = self.fluid.saturated_vapor_keyed_output(CP.iDmass)
else:
rho = self.rho[i]
self.mass_rate[i] = tp.gas_release_rate(
self.P[i],
self.p_back,
rho,
cpcv,
self.CD,
self.D_orifice**2 / 4 * math.pi,
)
elif input["valve"]["type"] == "hem_release":
if input["valve"]["flow"] == "filling":
raise ValueError(
"Filling flow not supported for HEM release valve type"
)
else:
self.mass_rate[i] = tp.hem_release_rate(
self.P[i],
self.p_back,
self.CD,
self.D_orifice**2 / 4 * math.pi,
self.fluid,
)
elif input["valve"]["type"] == "controlvalve":
Cv = tp.cv_vs_time(
self.Cv,
self.time_array[i],
self.valve_time_constant,
self.valve_characteristic,
)
if input["valve"]["flow"] == "filling":
Z = self.res_fluid.compressibility_factor()
MW = self.MW
k = self.res_fluid.cp0molar() / (self.res_fluid.cp0molar() - 8.314)
self.mass_rate[i] = -tp.control_valve(
self.p_back, self.P[i], self.T0, Z, MW, k, Cv
)
else:
Z = Z
MW = self.MW
# For NEM: use gas temperature for discharge
T_discharge = self.T_gas[i] if self.non_equilibrium and self.m_gas[i] > 1e-6 else self.T_fluid[i]
self.mass_rate[i] = tp.control_valve(
self.P[i], self.p_back, T_discharge, Z, MW, cpcv, Cv
)
elif input["valve"]["type"] == "psv":
# For NEM: use gas temperature for discharge
T_discharge = self.T_gas[i] if self.non_equilibrium and self.m_gas[i] > 1e-6 else self.T_fluid[i]
self.mass_rate[i] = tp.relief_valve(
self.P[i],
self.p_back,
self.Pset,
self.blowdown,
cpcv,
self.CD,
T_discharge,
Z,
self.MW,
self.D_orifice**2 / 4 * math.pi,
)
if (
"end_pressure" in self.input["valve"]
and self.input['valve']['flow']=='filling' and self.P[i] > self.input["valve"]["end_pressure"]
):
massflow_stop_switch = 1
elif (
"end_pressure" in self.input["valve"]
and self.input['valve']['flow']=='discharge' and self.P[i] < self.input["valve"]["end_pressure"]
):
massflow_stop_switch = 1
else:
massflow_stop_switch = 0
if massflow_stop_switch:
self.mass_rate[i] = 0
self.isrun = True
if input["valve"]["type"] == "relief":
idx_max = self.mass_rate.argmax()
# Smooth peak value by averaging with neighbors (avoid array index out of bounds)
if 0 < idx_max < len(self.mass_rate) - 1:
self.mass_rate[idx_max] = (
self.mass_rate[idx_max - 1] + self.mass_rate[idx_max + 1]
) / 2
self.relief_area[idx_max] = (
self.relief_area[idx_max - 1] + self.relief_area[idx_max + 1]
) / 2
# print("Relief area:", 2*math.sqrt(max(relief_area[1:])/math.pi), max(self.mass_rate))
[docs]
def get_dataframe(self):
"""
Export simulation results to pandas DataFrame for analysis and archiving.
Collects all time-series results from the simulation and organizes them
into a pandas DataFrame with labeled columns. Useful for exporting to
CSV/Excel, post-processing, or custom plotting.
Returns
-------
pd.DataFrame
DataFrame with simulation results. Columns include:
- Time (s)
- Pressure (bar)
- Fluid temperature (°C)
- Wall temperature (°C)
- Vent temperature (°C)
- Fluid enthalpy (J/kg)
- Fluid entropy (J/(kg·K))
- Fluid internal energy (J/kg)
- Discharge mass rate (kg/s)
- Fluid mass (kg)
- Fluid density (kg/m³)
- Inner heat transfer coefficient (W/(m²·K))
- Internal heat flux (W/m²)
- External heat flux (W/m²)
- Inner wall temperature (°C)
- Outer wall temperature (°C)
Notes
-----
Only returns data if simulation has been run (self.isrun == True).
Temperature values are converted from Kelvin to Celsius.
Pressure values are converted from Pa to bar.
"""
if self.isrun == True:
df = pd.DataFrame(self.time_array, columns=["Time (s)"])
df.insert(1, "Pressure (bar)", self.P / 1e5, True)
df.insert(2, "Fluid temperature (oC)", self.T_fluid - 273.15, True)
df.insert(3, "Wall temperature (oC)", self.T_vessel - 273.15, True)
df.insert(4, "Vent temperature (oC)", self.T_vent - 273.15, True)
df.insert(5, "Fluid enthalpy (J/kg)", self.H_mass, True)
df.insert(6, "Fluid entropy (J/kg K)", self.S_mass, True)
df.insert(7, "Fluid internal energy (J/kg)", self.U_mass, True)
df.insert(8, "Discharge mass rate (kg/s)", self.mass_rate, True)
df.insert(9, "Fluid mass (kg)", self.mass_fluid, True)
df.insert(10, "Fluid density (kg/m3)", self.rho, True)
df.insert(
11, "Inner heat transfer coefficient (W/m2 K)", self.h_inside, True
)
df.insert(
12,
"Internal heat flux (W/m2)",
self.Q_inner / self.surf_area_inner,
True,
)
df.insert(
13,
"External heat flux (W/m2)",
self.Q_outer / self.surf_area_outer,
True,
)
df.insert(
14, "Inner wall temperature (oC)", self.T_inner_wall - 273.15, True
)
df.insert(
13, "Outer wall temperature (oC)", self.T_outer_wall - 273.15, True
)
return df
[docs]
def plot(self, filename=None, verbose=True):
"""
Creating standard plots for the solved problem
Parameters
----------
filename : str
Saving plots to filename if provideed (optional)
verbose : bool
Plotting on screen if True (optional)
"""
import pylab as plt
if filename != None:
plt.figure(1, figsize=(12, 7), dpi=300)
else:
plt.figure(1, figsize=(8, 6))
plt.subplot(221)
# For NEM: plot gas and liquid temperatures separately
if self.non_equilibrium:
plt.plot(self.time_array, self.T_gas - 273.15, "r", label="Gas")
plt.plot(self.time_array, self.T_liquid - 273.15, "b", label="Liquid")
else:
plt.plot(self.time_array, self.T_fluid - 273.15, "b", label="Fluid")
if "thermal_conductivity" not in self.input["vessel"].keys():
plt.plot(
self.time_array, self.T_vessel - 273.15, "g", label="Vessel wall dry"
)
if self.liquid_level.any() != 0:
plt.plot(
self.time_array,
self.T_vessel_wetted - 273.15,
# marker="o",
color="darkorange",
label="Vessel wall wetted",
)
if "thermal_conductivity" in self.input["vessel"].keys():
if "liner_thermal_conductivity" in self.input["vessel"].keys():
plt.plot(
self.time_array,
self.T_bonded_wall - 273.15,
"g",
label="Liner/composite",
)
plt.plot(
self.time_array,
self.T_bonded_wall_wetted - 273.15,
color="darkorange",
label="Liner/composite wetted",
)
plt.plot(
self.time_array, self.T_inner_wall - 273.15, "g--", label="Inner wall"
)
plt.plot(
self.time_array, self.T_outer_wall - 273.15, "g-.", label="Outer wall"
)
plt.plot(
self.time_array,
self.T_inner_wall_wetted - 273.15,
color="darkorange",
linestyle="--",
label="Inner wall wetted",
)
plt.plot(
self.time_array,
self.T_outer_wall_wetted - 273.15,
color="darkorange",
linestyle="-.",
label="Outer wall wetted",
)
if self.input["valve"]["flow"] == "discharge":
plt.plot(self.time_array, self.T_vent - 273.15, "r", label="Vent")
if "validation" in self.input:
if "temperature" in self.input["validation"]:
temp = self.input["validation"]["temperature"]
if "gas_mean" in temp:
plt.plot(
np.asarray(temp["gas_mean"]["time"]),
np.asarray(temp["gas_mean"]["temp"]) - 273.15,
"b.",
label="Gas mean",
)
if "gas_high" in temp:
plt.plot(
np.asarray(temp["gas_high"]["time"]),
np.asarray(temp["gas_high"]["temp"]) - 273.15,
"b-.",
label="Gas high",
)
if "gas_low" in temp:
plt.plot(
np.asarray(temp["gas_low"]["time"]),
np.asarray(temp["gas_low"]["temp"]) - 273.15,
"b--",
label="Gas low",
)
if "wall_mean" in temp:
plt.plot(
np.asarray(temp["wall_mean"]["time"]),
np.asarray(temp["wall_mean"]["temp"]) - 273.15,
"go",
label="Wall mean",
)
if "wall_high" in temp:
plt.plot(
np.asarray(temp["wall_high"]["time"]),
np.asarray(temp["wall_high"]["temp"]) - 273.15,
"g-.",
label="Wall high",
)
if "wall_inner" in temp:
plt.plot(
np.asarray(temp["wall_inner"]["time"]),
np.asarray(temp["wall_inner"]["temp"]) - 273.15,
"g+",
label="Inner wall",
)
if "wall_low" in temp:
plt.plot(
np.asarray(temp["wall_low"]["time"]),
np.asarray(temp["wall_low"]["temp"]) - 273.15,
"g-.",
label="Wall high",
)
if "wall_outer" in temp:
plt.plot(
np.asarray(temp["wall_outer"]["time"]),
np.asarray(temp["wall_outer"]["temp"]) - 273.15,
"gx",
label="Outer wall",
)
plt.legend(loc="best")
plt.xlabel("Time (seconds)")
plt.ylabel(r"Temperature ($^\circ$C)")
plt.subplot(222)
plt.plot(self.time_array, self.P / 1e5, "b")
if "validation" in self.input:
if "pressure" in self.input["validation"]:
plt.plot(
np.asarray(self.input["validation"]["pressure"]["time"]),
self.input["validation"]["pressure"]["pres"],
"ko",
label="Experimental",
)
plt.legend(loc="best")
plt.xlabel("Time (seconds)")
plt.ylabel("Pressure (bar)")
plt.subplot(223)
plt.plot(
self.time_array, self.q_outer / 1000, "r", label="External heat flux (dry)"
)
plt.plot(
self.time_array,
self.q_inner / 1000,
color="darkorange",
label="Internal heat flux (dry)",
)
if self.liquid_level.any() != 0:
plt.plot(
self.time_array,
self.q_outer_wetted / 1000,
"r--",
label="External heat flux (wetted)",
)
plt.plot(
self.time_array,
self.q_inner_wetted / 1000,
color="darkorange",
linestyle="--",
label="Internal heat flux (wetted)",
)
plt.legend(loc="best")
plt.xlabel("Time (seconds)")
plt.ylabel("Heat flux (kW/m$^2$)")
plt.subplot(224)
plt.plot(self.time_array, self.mass_rate, "b", label="Mass flow (kg/s)")
plt.plot(self.time_array, self.liquid_level, "g", label="Liquid level (m)")
plt.legend(loc="best")
plt.xlabel("Time (seconds)")
plt.ylabel("Vent rate (kg/s) / Liquid level (m)")
plt.tight_layout()
if filename != None:
plt.savefig(filename + "_main.png")
if verbose:
plt.show()
return
[docs]
def plot_envelope(self, filename=None, verbose=True):
"""
Creating standard plots for the solved problem
Parameters
----------
filename : str
Saving plots to filename if provideed (optional)
verbose : bool
Plotting on screen if True (optional)
"""
import pylab as plt
if filename != None:
plt.figure(2, figsize=(12, 7), dpi=300)
else:
plt.figure(2, figsize=(8, 6))
self.fluid.build_phase_envelope("None")
PE = self.fluid.get_phase_envelope_data()
plt.plot(PE.T, PE.p, "-", label="HEOS Phase Envelope", color="g")
plt.plot(self.T_fluid, self.P, "-.", label="P/T fluid trajectory", color="b")
plt.plot(self.T_fluid[0], self.P[0], "o", label="Start", color="b")
plt.plot(self.T_fluid[-1], self.P[-1], ".", label="End", color="b")
plt.xlabel("Temperature [K]")
plt.ylabel("Pressure [Pa]")
plt.legend(loc="best")
plt.tight_layout()
if filename != None:
plt.savefig(filename + "_envelope.png")
if verbose:
plt.show()
[docs]
def plot_tprofile(self, filename=None, verbose=True):
"""
Creating standard plots for the solved problem
Parameters
----------
filename : str
Saving plots to filename if provideed (optional)
verbose : bool
Plotting on screen if True (optional)
"""
# Add some checks if the profile has been constructed
# return some
import pylab as plt
import numpy as np
if filename != None:
plt.figure(3, figsize=(8, 6))
else:
plt.figure(3, figsize=(8, 6))
X, Y = np.meshgrid(self.z * 1e3, self.time_array[:-1])
x0 = self.z[0] * 1e3
x1 = self.z[-1] * 1e3
y0 = self.time_array[0]
y1 = self.time_array[-1]
# plt.subplot(211)
plt.contourf(Y, X, np.asarray(self.temp_profile), origin="lower", levels=20)
# plt.imshow(np.asarray(self.temp_profile).T, aspect = 'auto', extent=(y0,y1,x0,x1), origin='lower')
plt.colorbar(label="Temperature (K)")
plt.xlabel("Time (s)")
plt.ylabel("z (mm)")
# add title with descriptive text for z axis
if filename != None:
plt.savefig(filename + "_tprofile1.png", dpi=300)
if filename != None:
plt.figure(4, figsize=(8, 6))
else:
plt.figure(4, figsize=(8, 6))
if verbose:
plt.show()
n = math.floor(len(self.time_array) / 15)
for i in range(len(self.time_array[::n])):
plt.plot(
self.temp_profile[::n][i],
self.z * 1e3,
label=f"t = {int(self.time_array[::n][i])} s.",
)
plt.legend(loc="best")
plt.ylabel("z (mm)")
plt.xlabel("Temperature (K)")
plt.title("Temperature distribution")
if filename != None:
plt.savefig(filename + "_tprofile2.png", dpi=300)
if verbose:
plt.show()
[docs]
def analyze_rupture(self, filename=None):
"""
Analyze vessel rupture potential under fire exposure conditions.
Performs a simplified rupture analysis by calculating vessel wall temperatures
under external fire heat load and comparing von Mises stress (from internal
pressure) against temperature-dependent allowable tensile stress (ATS).
Determines if and when vessel rupture may occur.
The analysis:
1. Calculates wall temperature evolution under fire exposure
2. Computes von Mises equivalent stress from internal pressure
3. Evaluates temperature-dependent material strength (ATS)
4. Identifies rupture time when stress exceeds strength
5. Generates diagnostic plots
Parameters
----------
filename : str, optional
Base filename for saving plots. If None, plots are displayed on screen.
Generates two plot files:
- {filename}_peak_wall_temp.png
- {filename}_ATS_vonmises.png
Returns
-------
None
Results are stored in instance variables:
- self.rupture_time : float or None
Time when rupture occurs [s], or None if no rupture predicted
- self.peak_times : ndarray
Time array for rupture analysis [s]
- self.von_mises : ndarray
von Mises equivalent stress history [Pa]
- self.ATS_wetted : ndarray
Allowable tensile stress for wetted region [Pa]
- self.ATS_unwetted : ndarray
Allowable tensile stress for unwetted region [Pa]
- self.peak_T_wetted : ndarray
Wall temperature history for wetted region [K]
- self.peak_T_unwetted : ndarray
Wall temperature history for unwetted region [K]
Notes
-----
Requires rupture analysis parameters in input:
- self.rupture_fire: Fire type (e.g., 'api_pool', 'scandpower_jet')
- self.rupture_material: Material type for ATS calculation
The method uses a simplified lumped capacitance model for wall heating.
Time step for rupture analysis is fixed at 10 seconds.
"""
from hyddown.materials import steel_Cp, ATS, von_mises
from hyddown import fire
pres = lambda x: np.interp(x, self.time_array, self.P)
q_unwetted = lambda x: np.interp(x, self.time_array, self.q_inner)
q_wetted = lambda x: np.interp(x, self.time_array, self.q_inner_wetted)
T0_unwetted = self.T_vessel[0]
T0_wetted = self.T_vessel_wetted[0]
thk = self.thickness
rho = self.vessel_density
inner_diameter = self.diameter
dt = 10
max_time = self.time_array[-1]
tsteps = int(max_time / dt)
T_wetted_wall = np.zeros(tsteps + 1)
T_unwetted_wall = np.zeros(tsteps + 1)
T_wetted_wall[0] = T0_wetted
T_unwetted_wall[0] = T0_unwetted
peak_times = np.zeros(tsteps + 1)
peak_times[0] = 0
for i in range(tsteps):
peak_times[i + 1] = peak_times[i] + dt
q_fire_wetted = fire.sb_fire(T_wetted_wall[i], self.rupture_fire)
q_fire_unwetted = fire.sb_fire(T_unwetted_wall[i], self.rupture_fire)
T_wetted_wall[i + 1] = T_wetted_wall[i] + (
q_fire_wetted - q_wetted(peak_times[i])
) * dt / (thk * rho * steel_Cp(T_wetted_wall[i], self.rupture_material))
T_unwetted_wall[i + 1] = T_unwetted_wall[i] + (
q_fire_unwetted - q_unwetted(peak_times[i])
) * dt / (thk * rho * steel_Cp(T_unwetted_wall[i], self.rupture_material))
ATS_wetted = np.array([ATS(T, self.rupture_material, k_s=self.rupture_k_s) for T in T_wetted_wall])
ATS_unwetted = np.array(
[ATS(T, self.rupture_material, k_s=self.rupture_k_s) for T in T_unwetted_wall]
)
von_mises_wetted = von_mises_unwetted = np.array(
[von_mises(pres(time), inner_diameter, thk) for time in peak_times]
)
self.peak_times = peak_times
self.von_mises = von_mises_unwetted
self.ATS_unwetted = ATS_unwetted
self.ATS_wetted = ATS_wetted
self.peak_T_wetted = T_wetted_wall
self.peak_T_unwetted = T_unwetted_wall
if np.all(ATS_unwetted > von_mises_unwetted) == True:
self.rupture_time = None
# print("No rupture")
elif np.all(ATS_unwetted < von_mises_unwetted) == True:
self.rupture_time = 0
# print("Rupture at time=0")
else:
rupture_idx = np.where(ATS_unwetted < von_mises_unwetted)[0][0]
self.rupture_time = (
peak_times[rupture_idx - 1] + peak_times[rupture_idx]
) / 2
# print("Rupture time +/- 5 s:", self.rupture_time)
# print("Rupture pressure (bar)", pres(peak_times[rupture_idx - 1]))
from matplotlib import pyplot as plt
plt.figure()
if np.any(self.liquid_level > 0):
plt.plot(peak_times, T_wetted_wall - 273.15, label="T wetted wall")
plt.plot(peak_times, T_unwetted_wall - 273.15, label="T unwetted wall")
plt.xlabel("Time (s)")
plt.ylabel("Wall temperature (C)")
plt.legend(loc="best")
if filename is not None:
plt.savefig(filename + "_peak_wall_temp.pdf")
plt.figure()
plt.plot(
peak_times,
np.array([pres(time) for time in peak_times]) / 1e5,
label="Pressure",
)
plt.xlabel("Time (s)")
plt.ylabel("Pressure (bar)")
plt.legend(loc="best")
if filename is not None:
plt.savefig(filename + "_peak_pressure.pdf")
plt.figure()
plt.plot(peak_times, von_mises_wetted / 1e6, label="von Mises stress")
plt.plot(peak_times, ATS_unwetted / 1e6, label="ATS unwetted wall")
if np.any(self.liquid_level > 0):
plt.plot(peak_times, ATS_wetted / 1e6, label="ATS wetted wall")
plt.xlabel("Time (s)")
plt.ylabel("Allowable Tensile Strength / von Mises Stress (MPa)")
plt.legend(loc="best")
if filename is not None:
plt.savefig(filename + "_ATS_vonmises.pdf")
if filename is None:
plt.show()
def __str__(self):
return "HydDown vessel filling/depressurization class"
[docs]
def generate_report(self):
"""
Generating a report summarising key features for the problem solved.
Can be used for e.g. case studies, problem optimisation (external) etc.
"""
report = {}
report["start_time"] = self.time_array[0]
report["end_time"] = self.time_array[-1]
# Pressure
report["max_pressure"] = max(self.P)
report["time_max_pressure"] = self.time_array[np.argmax(self.P)]
report["min_pressure"] = min(self.P)
report["time_min_pressure"] = self.time_array[np.argmin(self.P)]
# Temperatures
report["max_fluid_temp"] = max(self.T_fluid)
report["time_max_fluid_temp"] = self.time_array[np.argmax(self.T_fluid)]
report["min_fluid_temp"] = min(self.T_fluid)
report["time_min_fluid_temp"] = self.time_array[np.argmin(self.T_fluid)]
report["max_wall_temp"] = max(self.T_vessel)
report["time_max_wall_temp"] = self.time_array[np.argmax(self.T_vessel)]
report["min_wall_temp"] = min(self.T_vessel)
report["time_min_wall_temp"] = self.time_array[np.argmin(self.T_vessel)]
report["max_inner_wall_temp"] = max(self.T_inner_wall)
report["time_max_inner_wall_temp"] = self.time_array[
np.argmax(self.T_inner_wall)
]
report["min_inner_wall_temp"] = min(self.T_inner_wall)
report["time_min_inner_wall_temp"] = self.time_array[
np.argmin(self.T_inner_wall)
]
report["max_outer_wall_temp"] = max(self.T_outer_wall)
report["time_max_outer_wall_temp"] = self.time_array[
np.argmax(self.T_outer_wall)
]
report["min_outer_wall_temp"] = min(self.T_outer_wall)
report["time_min_outer_wall_temp"] = self.time_array[
np.argmin(self.T_outer_wall)
]
# Mass flows and inventory
report["max_mass_rate"] = max(self.mass_rate)
report["time_max_mass_rate"] = self.time_array[self.mass_rate.argmax()]
report["initial_mass"] = self.mass_fluid[0]
report["final_mass"] = self.mass_fluid[-1]
report["volume"] = self.vol
# Heat transfer (Q in W, q in W/m²)
# Track both max and min to capture extreme values in both directions
# (discharge: Q_inner > 0, filling: Q_inner < 0)
report["max_Q_inside"] = max(self.Q_inner)
report["time_max_Q_inside"] = self.time_array[np.argmax(self.Q_inner)]
report["min_Q_inside"] = min(self.Q_inner)
report["time_min_Q_inside"] = self.time_array[np.argmin(self.Q_inner)]
report["max_Q_outside"] = max(self.Q_outer)
report["time_max_Q_outside"] = self.time_array[np.argmax(self.Q_outer)]
report["min_Q_outside"] = min(self.Q_outer)
report["time_min_Q_outside"] = self.time_array[np.argmin(self.Q_outer)]
# Heat flux per unit area (use q arrays which store W/m²)
report["max_heat_flux_inside"] = max(self.q_inner)
report["time_max_heat_flux_inside"] = self.time_array[np.argmax(self.q_inner)]
report["min_heat_flux_inside"] = min(self.q_inner)
report["time_min_heat_flux_inside"] = self.time_array[np.argmin(self.q_inner)]
report["max_heat_flux_outside"] = max(self.q_outer)
report["time_max_heat_flux_outside"] = self.time_array[np.argmax(self.q_outer)]
report["min_heat_flux_outside"] = min(self.q_outer)
report["time_min_heat_flux_outside"] = self.time_array[np.argmin(self.q_outer)]
self.report = report