# 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)
# 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)
# 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)
[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 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
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
T_film = (
self.T_fluid[i - 1] + 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],
self.T_fluid[i - 1],
self.transport_fluid,
self.mass_rate[i - 1],
self.diameter,
)
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:
T_film = (
self.T_fluid[i - 1] + self.T_inner_wall[i - 1]
) / 2
try:
self.transport_fluid.update(
CP.PT_INPUTS, self.P[i - 1], T_film
)
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],
self.T_fluid[i - 1],
self.transport_fluid,
)
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
# Q = A * h * (T_wall - T_fluid)
self.Q_inner[i] = (
(self.surf_area_inner - wetted_area)
* hi
* (self.T_inner_wall[i - 1] - self.T_fluid[i - 1])
)
self.q_inner[i] = hi * (
self.T_inner_wall[i - 1] - self.T_fluid[i - 1]
)
# Heat transfer from wetted wall (liquid side) to fluid
# Uses different heat transfer coefficient (hiw) for liquid contact
self.Q_inner_wetted[i] = (
wetted_area
* hiw
* (self.T_inner_wall_wetted[i - 1] - self.T_fluid[i - 1])
)
# 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
)
if 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]
hi = tp.h_inner(
L,
self.T_fluid[i - 1],
T_wall_inner,
self.P[i - 1],
self.species,
)
self.h_inside[i] = hi
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
self.Q_inner[i] = self.scaling * (
(self.surf_area_inner - wetted_area)
* hi
* (T_wall_inner - self.T_fluid[i - 1])
)
self.q_inner[i] = hi * (T_wall_inner - self.T_fluid[i - 1])
self.Q_inner_wetted[i] = self.scaling * (
wetted_area * hiw * (T_wall_inner_wetted - self.T_fluid[i - 1])
)
self.q_inner_wetted[i] = hiw * (
T_wall_inner_wetted - self.T_fluid[i - 1]
)
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]
# Not pretty if-statement and a hack for fire relief area estimation. Most cases go directly to the first ...else... clause
if 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()
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
)
Z = self.fluid.saturated_vapor_keyed_output(CP.iZ)
else:
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:
if 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
self.mass_rate[i] = tp.control_valve(
self.P[i], self.p_back, self.T_fluid[i], Z, MW, cpcv, Cv
)
elif input["valve"]["type"] == "psv":
self.mass_rate[i] = tp.relief_valve(
self.P[i],
self.p_back,
self.Pset,
self.blowdown,
cpcv,
self.CD,
self.T_fluid[i],
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)
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) for T in T_wetted_wall])
ATS_unwetted = np.array(
[ATS(T, self.rupture_material) 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(1)
if self.liquid_level.all() != 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.png",
)
plt.figure(2)
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")
plt.figure(3)
plt.plot(peak_times, von_mises_wetted / 1e6, label="von Mises stress")
plt.plot(peak_times, ATS_unwetted / 1e6, label="ATS unwetted wall")
if self.liquid_level.all() != 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.png",
)
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