#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2024 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
#################################################################################
"""
IDAES 0D Fixed Bed Reactor model.
"""
# Import Pyomo libraries
from pyomo.environ import (
Constraint,
Var,
value,
check_optimal_termination,
units as pyunits,
)
from pyomo.dae import DerivativeVar
from pyomo.common.config import ConfigBlock, ConfigValue, In
# Import IDAES cores
from idaes.core import (
declare_process_block_class,
EnergyBalanceType,
UnitModelBlockData,
)
from idaes.core.util.config import is_physical_parameter_block
from idaes.core.util.exceptions import InitializationError
from idaes.core.util.constants import Constants as constants
import idaes.logger as idaeslog
from idaes.core.util import scaling as iscale
from idaes.core.solvers import get_solver
__author__ = "Chinedu Okoli, Andrew Lee"
# Assumptions:
# Only Vap and Sol phases, which are explicitly named as such
# Perfect mixing in Vap phase
# Static solid phase
[docs]
@declare_process_block_class("FixedBed0D")
class FixedBed0DData(UnitModelBlockData):
"""
0D Fixed Bed Reactor Model Class
"""
CONFIG = ConfigBlock()
CONFIG.declare(
"dynamic",
ConfigValue(
domain=In([True]),
default=True,
description="Dynamic model flag - must be True",
doc="""Indicates whether this model will be dynamic or not,
**default** = True. Fixed beds must be dynamic.""",
),
)
CONFIG.declare(
"has_holdup",
ConfigValue(
default=True,
domain=In([True]),
description="Holdup construction flag - must be True",
doc="""Indicates whether holdup terms should be constructed or not.
**default** - True. Fixed bed reactors must be dynamic, thus this must be
True.""",
),
)
CONFIG.declare(
"energy_balance_type",
ConfigValue(
default=EnergyBalanceType.enthalpyTotal,
domain=In(EnergyBalanceType),
description="Energy balance construction flag",
doc="""Indicates what type of energy balance should be constructed,
**default** - EnergyBalanceType.enthalpyTotal.
**Valid values:** {
**EnergyBalanceType.none** - exclude energy balances,
**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material,
**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase,
**EnergyBalanceType.energyTotal** - single energy balance for material,
**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
),
)
CONFIG.declare(
"gas_property_package",
ConfigValue(
default=None,
domain=is_physical_parameter_block,
description="Property package to use for gas phase",
doc="""Property parameter object used to define property calculations
for the gas phase, **default** - None.
**Valid values:** {
**useDefault** - use default package from parent model or flowsheet,
**PhysicalParameterObject** - a PhysicalParameterBlock object.}""",
),
)
CONFIG.declare(
"gas_property_package_args",
ConfigValue(
default=None,
domain=is_physical_parameter_block,
description="Arguments to use for constructing gas phase "
"property packages",
doc="""A ConfigBlock with arguments to be passed to a gas phase
property block(s) and used when constructing these,
**default** - None.
**Valid values:** {
see property package for documentation.}""",
),
)
CONFIG.declare(
"solid_property_package",
ConfigValue(
default=None,
domain=is_physical_parameter_block,
description="Property package to use for solid phase",
doc="""Property parameter object used to define property calculations
for the solid phase, **default** - None.
**Valid values:** {
**useDefault** - use default package from parent model or flowsheet,
**PhysicalParameterObject** - a PhysicalParameterBlock object.}""",
),
)
CONFIG.declare(
"solid_property_package_args",
ConfigValue(
default=None,
domain=is_physical_parameter_block,
description="Arguments to use for constructing solid phase "
"property packages",
doc="""A ConfigBlock with arguments to be passed to a solid phase
property block(s) and used when constructing these,
**default** - None.
**Valid values:** {
see property package for documentation.}""",
),
)
CONFIG.declare(
"reaction_package",
ConfigValue(
default=None,
# Removed domain argument due to limitations in ReactionBase for
# heterogeneous systems
description="Reaction package to use for unit",
doc="""Reaction parameter object used to define reaction calculations,
**default** - None.
**Valid values:** {
**None** - no reaction package,
**ReactionParameterBlock** - a ReactionParameterBlock object.}""",
),
)
CONFIG.declare(
"reaction_package_args",
ConfigBlock(
implicit=True,
description="Arguments to use for constructing reaction packages",
doc="""A ConfigBlock with arguments to be passed to a reaction block(s)
and used when constructing these,
**default** - None.
**Valid values:** {
see reaction package for documentation.}""",
),
)
[docs]
def build(self):
# Call UnitModel.build to setup dynamics
super().build()
# Get units meta data from property packages (only solid needed)
units_meta_solid = (
self.config.solid_property_package.get_metadata().get_derived_units
)
# Build Gas Phase StateBlock
# This block only needed so gas conc. to rxn block is calculated
# Defined state is set to True so that the "sum(mole_frac)=1" eqn in
# the gas state block is deactivated.
self.gas = self.config.gas_property_package.state_block_class(
self.flowsheet().time,
parameters=self.config.gas_property_package,
defined_state=True,
)
# Build Solid Phase StateBlock
# As there is no solid flow, only need a single state block
# Defined state is set to True so that the "sum(mass_frac)=1" eqn in
# the solid state block is deactivated. This is done here as there is
# currently no way to deactivate the constraint at the initial time
# for batch systems (i.e. no inlet or outlet ports).
# The "sum(mass_frac)=1 for all t neq 0" eqn is written in the
# unit model instead
self.solids = self.config.solid_property_package.state_block_class(
self.flowsheet().time,
parameters=self.config.solid_property_package,
defined_state=True,
)
tmp_dict = dict(**self.config.reaction_package_args)
tmp_dict["gas_state_block"] = self.gas
tmp_dict["solid_state_block"] = self.solids
tmp_dict["parameters"] = self.config.reaction_package
self.reactions = self.config.reaction_package.reaction_block_class(
self.flowsheet().time,
doc="Reaction properties in control volume",
**tmp_dict,
)
# Volume of reactor
self.bed_diameter = Var(
initialize=1, doc="Reactor diameter", units=units_meta_solid("length")
)
self.bed_height = Var(
initialize=1, doc="Bed length", units=units_meta_solid("length")
)
self.volume_bed = Var(
initialize=1.0,
doc="Volume of the reactor bed",
units=units_meta_solid("volume"),
)
@self.Constraint(doc="Calculating volume of the reactor bed")
def volume_bed_constraint(b):
return b.volume_bed == (
constants.pi * b.bed_height * (0.5 * b.bed_diameter) ** 2
)
# Solid phase material balance
# Volume of solid
self.volume_solid = Var(
self.flowsheet().config.time,
initialize=1.0,
doc="Solids phase volume including particle pores",
units=units_meta_solid("volume"),
)
@self.Constraint(
self.flowsheet().config.time, doc="Calculating solid phase volume"
)
def volume_solid_constraint(b, t):
return b.volume_solid[t] == (
b.volume_bed * (1 - b.solids[t].params.voidage)
)
# Accumulation equal to mass transfer/reaction
self.solids_material_holdup = Var(
self.flowsheet().time,
self.config.solid_property_package.component_list,
initialize=1,
doc="Solid phase component holdups",
units=units_meta_solid("mass"),
)
@self.Constraint(
self.flowsheet().config.time,
self.config.solid_property_package.component_list,
doc="Solid phase material holdup constraints",
)
def solids_material_holdup_constraints(b, t, j):
return b.solids_material_holdup[t, j] == (
b.volume_solid[t] * b.solids[t].get_material_density_terms("Sol", j)
)
self.solids_material_accumulation = DerivativeVar(
self.solids_material_holdup,
wrt=self.flowsheet().config.time,
doc="Solids material accumulation",
units=units_meta_solid("mass"),
)
# should be mass/time, to get around DAE limitations var is
# divided by time units when used
@self.Constraint(
self.flowsheet().config.time,
self.config.solid_property_package.component_list,
doc="Solid phase material accumulation constraints",
)
def solids_material_accumulation_constraints(b, t, j):
return b.solids_material_accumulation[t, j] / units_meta_solid(
"time"
) == b.volume_solid[t] * b.solids[t].params.mw_comp[j] * sum(
b.reactions[t].reaction_rate[r]
* b.config.reaction_package.rate_reaction_stoichiometry[r, "Sol", j]
for r in b.config.reaction_package.rate_reaction_idx
)
# Add solid mass variable and constraint for TGA tracking
self.mass_solids = Var(
self.flowsheet().config.time,
doc="Total mass of solids",
units=units_meta_solid("mass"),
)
@self.Constraint(
self.flowsheet().config.time, doc="Calculating total mass of solids"
)
def mass_solids_constraint(b, t):
return b.mass_solids[t] == b.volume_solid[t] * sum(
b.solids[t].get_material_density_terms("Sol", j)
for j in b.config.solid_property_package.component_list
)
# Sum of mass fractions at all time equals 1
@self.Constraint(
self.flowsheet().config.time, doc="Sum of mass fractions at all time"
)
def sum_component_constraint(b, t):
if t == b.flowsheet().config.time.first():
return Constraint.Skip
else:
return 1 == sum(
b.solids[t].mass_frac_comp[j]
for j in b.solids[t].params.component_list
)
if self.config.energy_balance_type != EnergyBalanceType.none:
# Solid phase energy balance
# Accumulation equal to heat transfer
self.solids_energy_holdup = Var(
self.flowsheet().time,
initialize=1,
doc="Solid phase energy holdup",
units=units_meta_solid("energy"),
)
@self.Constraint(
self.flowsheet().config.time,
doc="Solid phase energy holdup constraints",
)
def solids_energy_holdup_constraints(b, t):
return b.solids_energy_holdup[t] == (
self.volume_solid[t]
* pyunits.convert(
b.solids[t].get_energy_density_terms("Sol"),
to_units=units_meta_solid("energy")
/ units_meta_solid("volume"),
)
)
self.solids_energy_accumulation = DerivativeVar(
self.solids_energy_holdup,
wrt=self.flowsheet().config.time,
doc="Solids energy accumulation",
units=units_meta_solid("energy"),
)
# should be energy/time, to get around DAE limitations var is
# divided by time units when used
@self.Constraint(
self.flowsheet().config.time,
doc="Solid phase energy accumulation constraints",
)
def solids_energy_accumulation_constraints(b, t):
return b.solids_energy_accumulation[t] / units_meta_solid(
"time"
) == -sum(
b.reactions[t].reaction_rate[r]
* b.volume_solid[t]
* pyunits.convert(
b.reactions[t].dh_rxn[r],
to_units=units_meta_solid("energy_mole"),
)
for r in b.config.reaction_package.rate_reaction_idx
)
if self.config.energy_balance_type == EnergyBalanceType.none:
# Fix solids temperature to initial value for isothermal conditions
@self.Constraint(
self.flowsheet().config.time, doc="Isothermal solid phase constraint"
)
def isothermal_solid_phase(b, t):
if t == b.flowsheet().config.time.first():
return Constraint.Skip
else:
return b.solids[t].temperature == b.solids[0].temperature
[docs]
def initialize_build(
blk,
gas_phase_state_args=None,
solid_phase_state_args=None,
outlvl=idaeslog.NOTSET,
solver=None,
optarg=None,
):
"""
Initialization routine for FB0D unit.
Keyword Arguments:
gas_phase_state_args : a dict of arguments to be passed to the
property package(s) to provide an initial state for
initialization (see documentation of the specific
property package) (default = None).
solid_phase_state_args : a dict of arguments to be passed to the
property package(s) to provide an initial state for
initialization (see documentation of the specific
property package) (default = None).
outlvl : sets output level of initialization routine
optarg : solver options dictionary object (default=None, use
default solver options)
solver : str indicating which solver to use during
initialization (default = None, use default solver)
Returns:
None
"""
# Set solver options
init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
opt = get_solver(solver, optarg) # create solver
opt.options = optarg
# ---------------------------------------------------------------------
# Initialize control volume block
init_log.info("Initialize Thermophysical Properties")
# Initialize gas_phase block
gas_phase_flags = blk.gas.initialize(
outlvl=outlvl, optarg=optarg, solver=solver, state_args=gas_phase_state_args
)
# Initialize solid_phase properties block
solid_phase_flags = blk.solids.initialize(
outlvl=outlvl,
optarg=optarg,
solver=solver,
state_args=solid_phase_state_args,
)
init_log.info("Initialize reaction properties")
# Initialize reactions
blk.reactions.initialize(outlvl=outlvl, optarg=optarg, solver=solver)
init_log.info_high("Initialization Step 1 Complete.")
# ---------------------------------------------------------------------
# Solve unit
with idaeslog.solver_log(solve_log, outlvl) as slc:
res = opt.solve(blk, tee=slc.tee)
init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
# ---------------------------------------------------------------------
# Release Inlet state
blk.gas.release_state(gas_phase_flags, outlvl)
blk.solids.release_state(solid_phase_flags, outlvl)
if not check_optimal_termination(res):
raise InitializationError(
f"{blk.name} failed to initialize successfully. Please check "
f"the output logs for more information."
)
init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
def calculate_scaling_factors(self):
super().calculate_scaling_factors()
# scale some variables
if hasattr(self, "bed_diameter"):
if iscale.get_scaling_factor(self.bed_diameter) is None:
sf = 1 / value(self.bed_diameter)
iscale.set_scaling_factor(self.bed_diameter, 1e-1 * sf)
if hasattr(self, "bed_height"):
if iscale.get_scaling_factor(self.bed_height) is None:
sf = 1 / value(self.bed_height)
iscale.set_scaling_factor(self.bed_height, 1e-1 * sf)
if hasattr(self, "volume_bed"):
if iscale.get_scaling_factor(self.volume_bed) is None:
sf = 1 / value(
constants.pi * self.bed_height * (0.5 * self.bed_diameter) ** 2
)
iscale.set_scaling_factor(self.volume_bed, 1e-1 * sf)
# Scale some constraints
if hasattr(self, "volume_bed_constraint"):
for c in self.volume_bed_constraint.values():
iscale.constraint_scaling_transform(
c, iscale.get_scaling_factor(self.volume_bed), overwrite=False
)
if hasattr(self, "volume_solid_constraint"):
for t, c in self.volume_solid_constraint.items():
iscale.constraint_scaling_transform(
c, iscale.get_scaling_factor(self.volume_bed), overwrite=False
)
if hasattr(self, "solids_material_holdup_constraints"):
for (t, j), c in self.solids_material_holdup_constraints.items():
sf1 = iscale.get_scaling_factor(self.volume_bed)
sf2 = iscale.get_scaling_factor(self.solids[t].dens_mass_skeletal)
iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False)
if hasattr(self, "solids_material_accumulation_constraints"):
for (t, j), c in self.solids_material_accumulation_constraints.items():
iscale.constraint_scaling_transform(
c, iscale.get_scaling_factor(self.volume_bed), overwrite=False
)
if hasattr(self, "mass_solids_constraint"):
for t, c in self.mass_solids_constraint.items():
sf1 = iscale.get_scaling_factor(self.volume_bed)
sf2 = iscale.get_scaling_factor(self.solids[t].dens_mass_skeletal)
iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False)
if hasattr(self, "sum_component_constraint"):
# Get a single representative value in component list
for j in self.config.solid_property_package.component_list:
break
for t, c in self.sum_component_constraint.items():
iscale.constraint_scaling_transform(
c,
iscale.get_scaling_factor(self.solids[t].mass_frac_comp[j]),
overwrite=False,
)
if self.config.energy_balance_type != EnergyBalanceType.none:
if hasattr(self, "solids_energy_holdup_constraints"):
for t, c in self.solids_energy_holdup_constraints.items():
sf1 = iscale.get_scaling_factor(self.volume_bed)
sf2 = iscale.get_scaling_factor(self.solids[t].enth_mass)
iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False)
if hasattr(self, "solids_energy_accumulation_constraints"):
for t, c in self.solids_energy_accumulation_constraints.items():
sf1 = iscale.get_scaling_factor(self.volume_bed)
sf2 = iscale.get_scaling_factor(self.solids[t].enth_mass)
iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False)
if self.config.energy_balance_type == EnergyBalanceType.none:
if hasattr(self, "isothermal_solid_phase"):
for t, c in self.isothermal_solid_phase.items():
iscale.constraint_scaling_transform(
c,
iscale.get_scaling_factor(self.solids[t].temperature),
overwrite=False,
)
