Index

PackedBed

An implementation of the packed bed thermal energy storage model described by Battisti et al.¹.

Assumptions:

• Constant solid phase density
• Constant wall thermal properties
• Constant axial spacing
• Constant diameter
• Constant temperature exterior wall boundary condition

Warning

The model was intended for use with sCO₂ as the heat transfer fluid. The fluid attribute can be set to any CoolProp.AbstractState object, but the code has not been tested with any other fluid.

By default, sco2es.PackedBed uses CoolProp to calculate the properties of CO₂ and assumes that the solid particles are made of Alumina, for which empirical correlations are used to calculate relevant properties. To use a material other than alumina, create a class matching the protocol SolidProperties and set the solid attribute.

Info

Arrays are annotated with dimensions as follows:

• N: Number of time steps.
• Z: Number of axial nodes.
• W: Number of wall nodes.

---

1. F. Battisti, L. de Araujo Passos, and A. da Silva, “Performance mapping of packed-bed thermal energy storage systems for concentrating solar-powered plants using supercritical carbon dioxide,” Applied Thermal Engineering, vol. 183, p. 116032, 2021. ↩

Attributes:

Name	Type	Description
axial_nodes	int	Number of nodes in the axial direction.
wall_nodes	int	Number of nodes in the axial or radial directions for the walls and lids, respectively.
z	Annotated[NDArray[float], Z]	Axial positions of the packed bed nodes from the charging inlet in the direction of flow [m].
z_top_lid	Annotated[NDArray[float], W]	Axial positions of top lid nodes from charging inlet in direction of flow [m].
z_bottom_lid	Annotated[NDArray[float], W]	Axial positions of bottom lid nodes from charging inlet in direction of flow [m].
r_wall	Annotated[NDArray[float], W]	Radial positions of wall nodes from centerline [m].
T_f	Annotated[NDArray[float], (N, Z)]	Fluid temperature for bed nodes [K].
T_s	Annotated[NDArray[float], (N, Z)]	Solid temperature for bed nodes [K].
T_wall	Annotated[NDArray[float], (N, Z, W)]	Temperature of the wall nodes [K].
T_top_lid	Annotated[NDArray[float], (N, W)]	Temperature of the top lid nodes [K].
T_bottom_lid	Annotated[NDArray[float], (N, W)]	Temperature of the bottom lide nodes [K].
i_f	Annotated[NDArray[float], Z]	Fluid enthalpy at the node centers [J/kg].
cp_f	Annotated[NDArray[float], (N, Z)]	Fluid specific heat capacity at the node centers [J/kg⋅K].
k_f	Annotated[NDArray[float], Z]	Fluid thermal conductivity at the node centers [W/m⋅K].
rho_f	Annotated[NDArray[float], Z]	Fluid density at the node centers [kg/m3].
A_cs	float	Cross-sectional area of the packed bed [m2].
V_node	float	Volume of a packed bed node [m3].
r_bound	Annotated[NDArray[float], W + 1]	Radial positions of wall cell boundaries from centerline [m].
V_wall	Annotated[NDArray[float], W]	Volume of the wall cells [m3].
A_wall_z	Annotated[NDArray[float], W]	Surface area of the wall cell boundary in the axial direction [m2].
A_wall_r	Annotated[NDArray[float], W + 1]	Surface area of the wall cell boundary in the radial direction [m2].

Parameters:

Name	Type	Description	Default
T_initial	float	Discharge temperature [K].	required
P	float	Initial pressure [Pa].	required
L	float	Domain length [m].	required
D	float	Internal tank diameter [m].	required
d	float	Particle diameter [m].	required
eps	float	Void fraction.	required
T_env	float	Temperature of the environment [K].	required
t_wall	ArrayLike	Thickness of each wall layer [m].	required
k_wall	ArrayLike	Thermal conductivity of each wall layer [W/m⋅K].	required
rho_wall	ArrayLike	Density of each wall layer [kg/m3].	required
cp_wall	ArrayLike	Specific heat capacity of each wall layer [J/kg⋅K].	required
axial_nodes	int	Number of axial nodes.	100
wall_layer_nodes	int | ArrayLike	Number of nodes for each wall layer.	10
solid	SolidProperties | None	Solid properties interface.	None

fluid = CP.AbstractState('BICUBIC&HEOS', 'CO2') class-attribute instance-attribute

CoolProp.AbstractState object for accessing fluid properties. By default, CO₂ properties are tabulated for the Helmholtz-based equation of state (HEOS) and bicubic interpolation is used for performance.

solid: SolidProperties = Alumina class-attribute instance-attribute

[SolidPropsInterface][sco2es.SolidPropsInterface] for calculating temperature-dependent properties of the solid phase.

max_iter: int = 100 class-attribute instance-attribute

Maximum number of iterations for the loop.

atol_T_f: float = 0.05 class-attribute instance-attribute

Absolute tolerance for fluid temperature [ºC].

atol_P: float = 0.1 class-attribute instance-attribute

Absolute tolerance for pressure [Pa].

atol_T_s: float = 0.05 class-attribute instance-attribute

Absolute tolerance for solid temperature [ºC].

rtol_T_wall: float = 0.0005 class-attribute instance-attribute

Relative tolerance for wall and lid temperatures.

rtol_i_f: float = 0.0001 class-attribute instance-attribute

Relative tolerance for fluid enthalpy.

rtol_rho_f: float = 0.001 class-attribute instance-attribute

Relative tolerance for fluid density.

rtol_m_dot: float = 0.001 class-attribute instance-attribute

Relative tolerance for mass flow rate.

rtol_h: float = 0.001 class-attribute instance-attribute

Relative tolerance for volumetric and wall heat transfer coefficients.

E_loss_total property

Total energy loss at each timestep.

load_case classmethod

Loads a packed bed simulation case file in YAML format.

Parameters:

Name	Type	Description	Default
case_file	str	Filepath to the simulation case file.	required

time_index

A function for retrieving the index of the simulation time step with the closest elapsed time to the given time.

advance

Advances the simulation until the charge/discharge stopping criterion is satisfied.

Parameters:

Name	Type	Description	Default
T_inlet	float	Inlet fluid temperature [K].	required
P_inlet	float	Pressure at the inlet [Pa].	required
m_dot_inlet	float	Mass flow rate of the inlet stream [kg/s].	required
t_max	float	Maximum time advancement [s].	12 * 60 * 60
T_outlet_stop	float	Outlet fluid temperature [K] condition for stopping.	None
dt	float	Time step [s].	10
discharge	bool	Flag indicating that the system is discharging.	False

Returns:

Type	Description
	Time at which the charging/discharging stop criterion was satisfied.

Raises:

Type	Description
ModelAssumptionError	Raised if the Biot number exceeds the acceptable threshold (\(Bi > 0.1\)).
StopCriterionError	Raised if the stop criterion are not met within the maximum allowed charging/discharging time (t_max).

step

Calculates the state of the packed bed at the next time step using an iterative algorithm.

Parameters:

Name	Type	Description	Default
T_inlet	float	Temperature of the inlet stream [K].	required
P_inlet	float	Pressure at the inlet [Pa].	required
m_dot_inlet	float	Mass flow rate of the inlet stream [kg/s].	required
dt		Time step [s].	required
discharge	bool	Flag indicating that the system is discharging.	False

Returns:

Type	Description
	Number of iterations required for convergence.

solve_fluid_solid_bed staticmethod

Solves for the fluid enthalpy and solid temperature for the next time step.

Parameters:

Name	Type	Description	Default
i_f_0		Fluid enthalpy for the previous time step [J/kg].	required
i_inlet		Fluid enthalpy at the inlet [J/kg].	required
g		Fluid temperature-enthalpy coupling factor [K⋅kg/J].	required
e_s_0		Solid internal energy for the previous time step [J/kg].	required
alpha1		First linearized parameter for solid internal energy.	required
alpha2		Second linearized parameter for solid internal energy.	required
rho_f		Fluid density estimates for the next time step [kg/m3].	required
rho_0		Fluid densities for the previous time step [kg/m3].	required
rho_s		Density of the solid [kg/m3].	required
P_intf		Pressure estimate at the interfaces for the next time step [Pa].	required
P_intf_0		Pressure at the interfaces for the previous time step [Pa].	required
m_dot		Mass flow rate estimate at node interfaces [kg/s].	required
T_wall		Wall interface temperature [K].	required
T_lid_inlet		Temperature of the lid at the inlet [K].	required
T_lid_outlet		Temperature of the lid at the outlet [K].	required
k_eff		Effective thermal conductivity [W/m⋅K]	required
h_wall		Wall heat transfer coefficient [W/m2⋅K].	required
h_v		Volumetric heat transfer coefficient [W/m3⋅K].	required
A_wall		Surface area of the node wall boundary [m2].	required
A_cs		Cross-sectional area [m2].	required
V		Node volume [m3].	required
eps		Void fraction [-].	required
dz		Axial node spacing [m].	required
dt		Time step [s].	required

Returns:

Name	Type	Description
i_f		Fluid enthalpy for the next time step [J/kg].
T_s		Solid temperature for the next time step [K].

solve_lid_temperature staticmethod

Solves for the lid temperature profile for the next time step. The boundary conditions are taken as:

• Heat loss to environment by conduction at i=0
• Heat loss to fluid by convection at i=m-1
• Insulated at radial surfaces

The reverse keyword can be used to reverse the axial orientation of the boundary conditions.

Parameters:

Name	Type	Description	Default
T_lid		Lid temperature for the previous time step [K].	required
T_f		Fluid temperature estimate for the next time step at convection boundary [K].	required
T_env		Environment temperature [K].	required
h_wall		Wall heat transfer coefficient [W/m2⋅K].	required
k		Thermal conductivity [W/m⋅K].	required
rho		Density [kg/m3].	required
cp		Specific heat capacity [J/kg⋅K].	required
z		Axial positions of node centers [m].	required
V		Node volumes [m3].	required
A		Area of lid [m2].	required
dt		Time step [s].	required
reverse		Flag for reversing boundary conditions.	False

solve_wall_temperature staticmethod

Solves for the lid temperature profile for the next time step. The boundary conditions are taken as:

• Heat loss to fluid by convection at j=0
• Heat loss to environment by conduction at j=m-1
• Insulated axially at i=0 and i=n-1

where the domain is discretized with n nodes in the axial direction and m in the radial direction, and the wall temperature T_wall has the shape (n, m).

Parameters:

Name	Type	Description	Default
T_wall		Wall temperature for the previous time step [K].	required
T_f		Fluid temperature estimate for the next time step [K].	required
T_env		Environment temperature [K].	required
h_wall		Wall heat transfer coefficient [W/m2⋅K].	required
k		Thermal conductivity [W/m K].	required
rho		Density [kg/m3].	required
cp		Specific heat capacity [J/kg⋅K].	required
r		Radial positions of node centers [m].	required
dr		Radial node width [m].	required
dz		Axial node spacing [m].	required
V		Node volumes [m3].	required
A_r		Area of node interfaces in the radial direction [m2].	required
A_z		Area of node interfaces in the axial direction [m2].	required
dt		Time step [s].	required

calculate_heat_transfer_coeffs

Calculates the relevant heat transfer coefficients for the packed bed.

Parameters:

Name	Type	Description	Default
m_dot		Mass flow rate [kg/s].	required
T_f		Temperature of the fluid [K].	required
k_f		Thermal conductivity of the fluid [W/m⋅K].	required
cp_f		Specific heat capacity of the fluid [J/kg⋅K].	required
mu_f		Dynamic viscosity of the fluid [Pa⋅s].	required
k_s		Thermal conductivity of the solid [W/m⋅K].	required
E_s		Emissivity of the solid.	required

Returns:

Name	Type	Description
k_eff		Effective thermal conductivity [W/m⋅K].
h_wall		Wall heat transfer coefficient [W/m2⋅K].
h_v		Volumetric heat transfer coefficient [W/m3⋅K].

calculate_fluid_enthalpy

Returns the fluid mass-specific enthalpy at each node.

Parameters:

Name	Type	Description	Default
T_f		Temperature [K].	required
P		Pressure [Pa].	required

calculate_fluid_props

Returns the fluid temperature, thermal conductivity, density, viscosity, and specific heat capacity at each node.

Parameters:

Name	Type	Description	Default
i_f		Mass-specific enthalpy of the fluid [J/kg].	required
P		Pressure [Pa].	required

Returns:

Name	Type	Description
T_f		Temperature of the fluid [K].
k_f		Thermal conductivity of the fluid [W/m⋅K].
rho_f		Density of the fluid [kg/m3].
mu_f		Dynamic viscosity of the fluid [Pa⋅s].
cp_f		Specific heat capacity of the fluid [J/kg⋅K].

biot_number staticmethod

Calculates the biot number (\(Bi\)) for a packed bed¹

\[ Bi = \frac{h_v d^2}{36(1 - \varepsilon) k_s} \]

---

1. F. Battisti, L. de Araujo Passos, and A. da Silva, “Performance mapping of packed-bed thermal energy storage systems for concentrating solar-powered plants using supercritical carbon dioxide,” Applied Thermal Engineering, vol. 183, p. 116032, 2021. ↩

Parameters:

Name	Type	Description	Default
h_v		Volumetric heat transfer coefficient, \(h_v\) [W/m3⋅K].	required
d		Particle diameter, \(d\) [m].	required
eps		Void fraction, \(\varepsilon\).	required
k_s		Thermal conductivity of the solid, \(k_s\) [W/m⋅K].	required

effective_film_thickness_ratio staticmethod

Calculates the ratio between the effective thickness of the fluid film adjacent to the surface of two solid particles and the particle diameter according to the interpolation of Kunii and Smith¹ for \(0.260 \leq \varepsilon \leq 0.476\)

\[ \phi = \phi_2 + (\phi_1 - \phi_2) \frac{\varepsilon - 0.260}{0.216} \]

where

\[ \phi_i = \frac{1}{2} \frac{\left(\frac{\kappa-1}{\kappa}\right)^2 \sin^2{\theta_i}} {\ln\left[\kappa - (\kappa - 1)\cos{\theta_i}\right] - \left(\frac{\kappa - 1}{\kappa}\right) (1 - \cos{\theta_i)}} - \frac{2}{3\kappa} \]

where \(\kappa = k_s/k_f\) and the interpolation bounds are given by \(\sin^2 \theta_i = \frac{1}{n_i}\) for \(n_1 = 1.5\) and \(n_2 = 4\sqrt{3}\).

---

1. D. Kunii and J. Smith, “Heat transfer characteristics of porous rocks,” AIChE Journal, vol. 6, no. 1, pp. 71–78, 1960. ↩

Parameters:

Name	Type	Description	Default
k_s		Thermal conductivity of the solid, \(k_s\) [W/m⋅K].	required
k_f		Thermal conductivity of the fluid, \(k_f\) [W/m⋅K].	required
eps		Void fraction, \(\varepsilon\).	required

void_radiative_heat_transfer_coeff staticmethod

Calculates the void-to-void radiative heat transfer coefficient for a packed bed given by the correlation of Yagi and Kunii¹

\[ h_{rv} = \frac{0.1952}{1 + \frac{\varepsilon (1-E_s)}{2E_s (1-\varepsilon)}} \left(\frac{T}{100}\right)^3 \]

---

1. S. Yagi and D. Kunii, “Studies on effective thermal conductivities in packed beds,” AIChE Journal, vol. 3, no. 3, pp. 373–381, 1957. ↩

Parameters:

Name	Type	Description	Default
T		Temperature, \(T\) [K].	required
eps		Void fraction, \(\varepsilon\).	required
E_s		Emissivity of the solid, \(E_s\).	required

surface_radiative_heat_transfer_coeff staticmethod

Calculates the surface-to-surface radiative heat transfer coefficient for a packed bed given by the correlation of Yagi and Kunii¹

\[ h_{rs} = \frac{0.1952E_s}{2-E_s} \left( \frac{T}{100}\right)^3 \]

---

1. S. Yagi and D. Kunii, “Studies on effective thermal conductivities in packed beds,” AIChE Journal, vol. 3, no. 3, pp. 373–381, 1957. ↩

Parameters:

Name	Type	Description	Default
T		Temperature, \(T\) [K].	required
E_s		Emissivity of the solid, \(E_s\).	required

effective_thermal_conductivity staticmethod

Approximates the effective thermal conductivity in a packed bed using the model of Kunii and Smith¹

\[ k_{eff} = k_f \left[\varepsilon\left(1+\beta \frac{h_{r\nu}d}{k_f}\right) + \frac{\beta (1-\varepsilon}{\frac{1}{\frac{1}{\phi}+\frac{h_{rs}d}{k_f}} +\gamma\left(\frac{k_f}{k_s}\right)} \right] \]

---

1. D. Kunii and J. Smith, “Heat transfer characteristics of porous rocks,” AIChE Journal, vol. 6, no. 1, pp. 71–78, 1960. ↩

Parameters:

Name	Type	Description	Default
k_f		Thermal conductivity of the fluid, \(k_f\) [W/m⋅K].	required
k_s		Thermal conductivity of the solid, \(k_s\) [W/m⋅K].	required
eps		Void fraction. \(\varepsilon\).	required
h_rv		Void-to-void radiative heat transfer coefficient, \(h_{rv}\) [W/m2⋅K].	required
h_rs		Surface-to-surface radiative heat transfer coefficient, \(h_{rs}\) [W/m2⋅K].	required
phi		Effective film thickness ratio, \(\phi\).	required
d		Solid particle diameter, \(d\) [m].	required
beta		Ratio of the average length between the centers of two neighboring solids to the mean particle diameter, \(\beta = 0.9\) (default).	0.9

volumetric_convective_heat_transfer_coeff staticmethod

Calculates the volumetric convective heat transfer coefficient as the product of the convective heat transfer coefficient for a spherical particle and the ratio between the particles total heat transfer area and the total volume

\[ h_v = h_{part} \frac{6(1-\varepsilon)}{d} \]

The particle convective heat transfer coefficient is given by Pfeffer's correlation¹

\[ h_{part} = 1.26 \left[\frac{1-(1-\varepsilon)^{5/3}}{W}\right]^{1/3} (c_{p_f} G)^{1/3} \left(\frac{k_f}{d}\right)^{2/3} \]

where \(W\) is given by

\[ W = 2-3(1-\varepsilon)^{1/3}+3(1-\varepsilon)^{5/3}-2(1-\varepsilon)^2 \]

and \(G\) is the effective mass flow rate per unit area given by

\[ G = \frac{4\dot{m}}{\varepsilon \pi D^2} \]

The lower limit of analytical heat transfer is given by

\[ h_{part} = \frac{2k_f}{d} \]

for a heated isothermal sphere in a quiescent fluid medium; therefore, the maximum of the two values is taken.

---

1. R. Pfeffer, “Heat and mass transport in multiparticle systems,” Industrial & Engineering Chemistry Fundamentals, vol. 3, no. 4, pp. 380–383, 1964. ↩

Parameters:

Name	Type	Description	Default
m_dot		Mass flow rate [kg/s].	required
k_f		Thermal conductivity of the fluid, \(k_f\) [W/m⋅K].	required
cp_f		Specific heat capacity of the fluid, \(c_{p_f}\) [J/kg⋅K].	required
eps		Void fraction, \(\varepsilon\).	required
d		Particle diameter, \(d\) [m].	required
D		Diameter [m].	required

conv_wall_heat_transfer_coeff staticmethod

Returns the convective heat transfer coefficient between the fluid and the wall

\[ h_{wall}^{cv} = \left( 2.58 Re_d^{1/3} Pr^{1/3} + 0.094 Re_d^{0.8} Pr^{0.4} \right) \frac{k_f}{d} \]

according to the correlation of Beek¹. The Reynolds number \(Re_d\) is defined as

\[ Re_d = \frac{\rho_f u_0 d}{\mu_f} \]

where \(u_0\) is the superficial velocity, or the velocity if no particles were present, given by

\[ u_0 = \frac{\dot{m}}{\rho_f A} \]

The Prandtl number \(Pr\) is defined as

\[ Pr = \frac{{c_p}_f \mu_f}{k_f} \]

---

1. J. Beek, “Design of packed catalytic reactors,” in Advances in Chemical Engineering, Elsevier, 1962, pp. 203–271. ↩

Parameters:

Name	Type	Description	Default
m_dot		Mass flow rate, \(\dot{m}\) [kg/s].	required
k_f		Thermal conductivity of the fluid, \(k_f\) [W/m⋅K].	required
cp_f		Specific heat capacity of the fluid, \({c_p}_f\) [J/kg⋅K].	required
mu_f		Dynamic viscosity of the fluid, \(\mu_f\) [Pa⋅s].	required
d		Particle diameter, \(d\) [m].	required
D		Diameter, \(D\) [m].	required

cond_rad_wall_heat_transfer_coeff staticmethod

Calculates the conductive and radiative wall heat transfer coefficient according to Ofuchi and Kunii¹

\[ h_{wall}^{cd,ra} = \frac{k_{eff}^{stag} k_{wall}^{stag}}{k_{eff}^{stag} - \frac{k_{wall}^{stag}}{2}} \]

where

\[ k_{eff}^{stag} = k_f \left[\varepsilon \left(1 + \frac{h_{rv}d}{k_f}\right) + \frac{1-\varepsilon}{\left(\frac{1}{\phi}+\frac{h_{rs}d}{k_f}\right)^{-1} + \frac{2}{3\kappa}}\right] \]

and

\[ k_{wall}^{stag} = k_f \left[\varepsilon_{wall} \left(2 + \frac{h_{rv}d}{k_f}\right) + \frac{1-\varepsilon_{wall}}{\left(\frac{1}{\phi_{wall}}+\frac{h_{rs}d}{k_f}\right)^{-1} + \frac{1}{3\kappa}}\right] \]

where the wall porosity, \(\varepsilon_{wall}\), is assumed to be \(0.4\) and \(\phi_{wall}\) is given by

\[ \phi_{wall} = \frac{1}{4} \frac{\left(\frac{\kappa-1}{\kappa}\right)^2}{\ln{\kappa} - \frac{\kappa-1}{\kappa}} - \frac{1}{3\kappa} \]

where \(\kappa = k_s/k_f\).

---

1. K. Ofuchi and D. Kunii, “Heat-transfer characteristics of packed beds with stagnant fluids,” International Journal of Heat and Mass Transfer, vol. 8, no. 5, pp. 749–757, 1965. ↩

Parameters:

Name	Type	Description	Default
k_f		Thermal conductivity of the fluid, \(k_f\) [W/m⋅K].	required
k_s		Thermal conductivity of the solid, \(k_s\) [W/m⋅K].	required
h_rv		void-to-void radiative heat transfer coefficient, \(h_{rv}\) [W/m2⋅K].	required
h_rs		Solid-to-solid radiative heat transfer coefficient, \(h_{rs}\) [W/m2⋅K].	required
eps		Void fraction, \(\varepsilon\).	required
d		Particle diameter, \(d\) [m].	required
phi		Effective film thickness ratio, \(\phi\).	required

pressure_drop staticmethod

Calculates the pressure drop using the modified Ergun's equation¹

\[ \Delta P = \frac{\Delta z G^2}{\rho_f d} \left[ \xi_1 \frac{(1-\epsilon)^2}{\epsilon^3 \psi^2} \frac{\mu_f}{Gd} + \xi_2 \frac{(1-\epsilon)}{\epsilon^3 \psi}\right] \]

---

1. I. Macdonald, M. El-Sayed, K. Mow, and F. Dullien, “Flow through porous media-the Ergun equation revisited,” Industrial & Engineering Chemistry Fundamentals, vol. 18, no. 3, pp. 199–208, 1979. ↩

Parameters:

Name	Type	Description	Default
dz	float	Axial position step size, \(\Delta z\) [m].	required
rho_f	float	Density of the fluid, \(\rho_f\) [kg/m3].	required
mu_f	float	Dynamic viscosity of the fluid, \(\mu_f\) [Pa⋅s].	required
G	float	Effective mass flow rate per unit cross-section, \(G\) [kg/m2⋅s].	required
eps	float	Void fraction, \(\varepsilon\).	required
d	float	Particle diameter, \(d\) [m].	required
psi	float	Sphericity, \(\psi\).	0.9
xi1	float	Viscous loss coefficient, \(\xi_1\).	180
xi2	float	Inertial loss coefficient, \(\xi_2\).	1.8