Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding support for inherent reactions in Mixer and Separator #1289

Merged
merged 5 commits into from
Nov 20, 2023
Merged

Adding support for inherent reactions in Mixer and Separator #1289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f187fff
Removing unnecessary 0 terms in material balances
Nov 15, 2023
ff5ac46
Testing split with inherent rxns
Nov 16, 2023
f24f0f1
Adding support for inherent rxns to mixer
Nov 16, 2023
9db90a6
Fixing pylint issues
Nov 16, 2023
4d428e5
Merge branch 'main' into mix_split_inherent_rxns
Nov 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 79 additions & 31 deletions idaes/models/unit_models/mixer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
from pyomo.environ import (
Block,
check_optimal_termination,
Constraint,
Param,
PositiveReals,
Reals,
Expand Down Expand Up @@ -635,6 +636,50 @@
# Let this pass for now with no units
flow_units = None

if mixed_block.include_inherent_reactions:
if mb_type == MaterialBalanceType.total:
raise ConfigurationError(
"Cannot do total flow mixing with inherent reaction; "
"problem is under-constrained. Please use a different "
"mixing type."
)

# Add extents of reaction and stoichiometric constraints
self.inherent_reaction_extent = Var(
self.flowsheet().time,
mixed_block.params.inherent_reaction_idx,
domain=Reals,
initialize=0.0,
doc="Extent of inherent reactions in outlet",
units=flow_units,
)

self.inherent_reaction_generation = Var(
self.flowsheet().time,
pc_set,
domain=Reals,
initialize=0.0,
doc="Generation due to inherent reactions in outlet",
units=flow_units,
)

@self.Constraint(
self.flowsheet().time,
pc_set,
)
def inherent_reaction_constraint(b, t, p, j):
if (p, j) in pc_set:
return b.inherent_reaction_generation[t, p, j] == (
sum(
mixed_block[t].params.inherent_reaction_stoichiometry[
r, p, j
]
* self.inherent_reaction_extent[t, r]
for r in mixed_block[t].params.inherent_reaction_idx
)
)
return Constraint.Skip

Check warning on line 681 in idaes/models/unit_models/mixer.py

View check run for this annotation

Codecov / codecov/patch

idaes/models/unit_models/mixer.py#L681 

Added line #L681 was not covered by tests

if mb_type == MaterialBalanceType.componentPhase:
# Create equilibrium generation term and constraints if required
if self.config.has_phase_equilibrium is True:
Expand All @@ -653,43 +698,35 @@
"thus does not support phase equilibrium.".format(self.name)
)

# Define terms to use in mixing equation
def phase_equilibrium_term(b, t, p, j):
if self.config.has_phase_equilibrium:
sd = {}
for r in pp.phase_equilibrium_idx:
if pp.phase_equilibrium_list[r][0] == j:
if pp.phase_equilibrium_list[r][1][0] == p:
sd[r] = 1
elif pp.phase_equilibrium_list[r][1][1] == p:
sd[r] = -1
else:
sd[r] = 0
else:
sd[r] = 0

return sum(
b.phase_equilibrium_generation[t, r] * sd[r]
for r in pp.phase_equilibrium_idx
)
else:
return 0

# Write phase-component balances
@self.Constraint(
self.flowsheet().time,
pc_set,
doc="Material mixing equations",
)
def material_mixing_equations(b, t, p, j):
return 0 == (
sum(
inlet_blocks[i][t].get_material_flow_terms(p, j)
for i in range(len(inlet_blocks))
rhs = sum(
inlet_blocks[i][t].get_material_flow_terms(p, j)
for i in range(len(inlet_blocks))
) - mixed_block[t].get_material_flow_terms(p, j)

if self.config.has_phase_equilibrium:
rhs += sum(
b.phase_equilibrium_generation[t, r]
for r in pp.phase_equilibrium_idx
if pp.phase_equilibrium_list[r][0] == j
and pp.phase_equilibrium_list[r][1][0] == p
) - sum(
b.phase_equilibrium_generation[t, r]
for r in pp.phase_equilibrium_idx
if pp.phase_equilibrium_list[r][0] == j
and pp.phase_equilibrium_list[r][1][1] == p
)
- mixed_block[t].get_material_flow_terms(p, j)
+ phase_equilibrium_term(b, t, p, j)
)

if mixed_block.include_inherent_reactions:
rhs += b.inherent_reaction_generation[t, p, j]

return 0 == rhs

elif mb_type == MaterialBalanceType.componentTotal:
# Write phase-component balances
Expand All @@ -699,7 +736,7 @@
doc="Material mixing equations",
)
def material_mixing_equations(b, t, j):
return 0 == sum(
rhs = sum(
sum(
inlet_blocks[i][t].get_material_flow_terms(p, j)
for i in range(len(inlet_blocks))
Expand All @@ -709,11 +746,20 @@
if (p, j) in pc_set
)

if mixed_block.include_inherent_reactions:
rhs += sum(
b.inherent_reaction_generation[t, p, j]
for p in mixed_block.phase_list
if (p, j) in pc_set
)

return 0 == rhs

elif mb_type == MaterialBalanceType.total:
# Write phase-component balances
@self.Constraint(self.flowsheet().time, doc="Material mixing equations")
def material_mixing_equations(b, t):
return 0 == sum(
rhs = sum(
sum(
sum(
inlet_blocks[i][t].get_material_flow_terms(p, j)
Expand All @@ -726,6 +772,8 @@
for p in mixed_block.phase_list
)

return 0 == rhs

elif mb_type == MaterialBalanceType.elementTotal:
raise ConfigurationError(
"{} Mixers do not support elemental "
Expand Down
153 changes: 107 additions & 46 deletions idaes/models/unit_models/separator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -932,21 +932,65 @@
t_ref = self.flowsheet().time.first()
mb_type = mixed_block[t_ref].default_material_balance_type()

# Get units from property package
units_meta = mixed_block.params.get_metadata()
flow_basis = mixed_block[
self.flowsheet().time.first()
].get_material_flow_basis()
if flow_basis == MaterialFlowBasis.molar:
flow_units = units_meta.get_derived_units("flow_mole")
elif flow_basis == MaterialFlowBasis.mass:
flow_units = units_meta.get_derived_units("flow_mass")
else:
# Let this pass for now with no units
flow_units = None

Check warning on line 946 in idaes/models/unit_models/separator.py

View check run for this annotation

Codecov / codecov/patch

idaes/models/unit_models/separator.py#L946 

Added line #L946 was not covered by tests

if mixed_block.include_inherent_reactions:
# Add extents of reaction and stoichiometric constraints
# TODO: It would be nice if there was a way to recognise that
# a total flow split does not require calculation of equilibrium
# (assuming T&P are constant) and turn of inherent reactions in
# that case, but that feature does not exist in property packages
self.inherent_reaction_extent = Var(
self.flowsheet().time,
self.outlet_idx,
mixed_block.params.inherent_reaction_idx,
domain=Reals,
initialize=0.0,
doc="Extent of inherent reactions in outlets",
units=flow_units,
)

self.inherent_reaction_generation = Var(
self.flowsheet().time,
self.outlet_idx,
pc_set,
domain=Reals,
initialize=0.0,
doc="Generation due to inherent reactions in outlets",
units=flow_units,
)

@self.Constraint(
self.flowsheet().time,
self.outlet_idx,
pc_set,
)
def inherent_reaction_constraint(b, t, o, p, j):
if (p, j) in pc_set:
return b.inherent_reaction_generation[t, o, p, j] == (
sum(
mixed_block[t].params.inherent_reaction_stoichiometry[
r, p, j
]
* self.inherent_reaction_extent[t, o, r]
for r in mixed_block[t].params.inherent_reaction_idx
)
)
return Constraint.Skip

Check warning on line 990 in idaes/models/unit_models/separator.py

View check run for this annotation

Codecov / codecov/patch

idaes/models/unit_models/separator.py#L990 

Added line #L990 was not covered by tests

if mb_type == MaterialBalanceType.componentPhase:
if self.config.has_phase_equilibrium is True:
# Get units from property package
units_meta = mixed_block.params.get_metadata()
flow_basis = mixed_block[
self.flowsheet().time.first()
].get_material_flow_basis()
if flow_basis == MaterialFlowBasis.molar:
flow_units = units_meta.get_derived_units("flow_mole")
elif flow_basis == MaterialFlowBasis.mass:
flow_units = units_meta.get_derived_units("flow_mass")
else:
# Let this pass for now with no units
flow_units = None

try:
self.phase_equilibrium_generation = Var(
self.flowsheet().time,
Expand All @@ -963,29 +1007,6 @@
"does not support phase equilibrium.".format(self.name)
)

# Define terms to use in mixing equation
def phase_equilibrium_term(b, t, o, p, j):
if self.config.has_phase_equilibrium:
sd = {}
sblock = mixed_block[t]
for r in b.config.property_package.phase_equilibrium_idx:
if sblock.params.phase_equilibrium_list[r][0] == j:
if sblock.params.phase_equilibrium_list[r][1][0] == p:
sd[r] = 1
elif sblock.params.phase_equilibrium_list[r][1][1] == p:
sd[r] = -1
else:
sd[r] = 0
else:
sd[r] = 0

return sum(
b.phase_equilibrium_generation[t, o, r] * sd[r]
for r in b.config.property_package.phase_equilibrium_idx
)
else:
return 0

@self.Constraint(
self.flowsheet().time,
self.outlet_idx,
Expand All @@ -994,12 +1015,27 @@
)
def material_splitting_eqn(b, t, o, p, j):
o_block = getattr(self, o + "_state")
return sf(t, o, p, j) * mixed_block[t].get_material_flow_terms(
p, j
) == (
o_block[t].get_material_flow_terms(p, j)
- phase_equilibrium_term(b, t, o, p, j)
)
lhs = sf(t, o, p, j) * mixed_block[t].get_material_flow_terms(p, j)

rhs = o_block[t].get_material_flow_terms(p, j)

if self.config.has_phase_equilibrium:
rhs += -sum(
b.phase_equilibrium_generation[t, o, r]
for r in b.config.property_package.phase_equilibrium_idx
if mixed_block[t].params.phase_equilibrium_list[r][0] == j
and mixed_block[t].params.phase_equilibrium_list[r][1][0] == p
) + sum(
b.phase_equilibrium_generation[t, o, r]
for r in b.config.property_package.phase_equilibrium_idx
if mixed_block[t].params.phase_equilibrium_list[r][0] == j
and mixed_block[t].params.phase_equilibrium_list[r][1][1] == p
)

if mixed_block.include_inherent_reactions:
rhs += b.inherent_reaction_generation[t, o, p, j]

return lhs == rhs

elif mb_type == MaterialBalanceType.componentTotal:

Expand All @@ -1011,16 +1047,27 @@
)
def material_splitting_eqn(b, t, o, j):
o_block = getattr(self, o + "_state")
return sum(

lhs = sum(
sf(t, o, p, j) * mixed_block[t].get_material_flow_terms(p, j)
for p in mixed_block.phase_list
if (p, j) in pc_set
) == sum(
)
rhs = sum(
o_block[t].get_material_flow_terms(p, j)
for p in o_block.phase_list
if (p, j) in pc_set
)

if mixed_block.include_inherent_reactions:
rhs += sum(
b.inherent_reaction_generation[t, o, p, j]
for p in o_block.phase_list
if (p, j) in pc_set
)

return lhs == rhs

elif mb_type == MaterialBalanceType.total:

@self.Constraint(
Expand All @@ -1030,14 +1077,16 @@
)
def material_splitting_eqn(b, t, o):
o_block = getattr(self, o + "_state")
return sum(

lhs = sum(
sum(
sf(t, o, p, j) * mixed_block[t].get_material_flow_terms(p, j)
for j in mixed_block.component_list
if (p, j) in pc_set
)
for p in mixed_block.phase_list
) == sum(
)
rhs = sum(
sum(
o_block[t].get_material_flow_terms(p, j)
for j in mixed_block.component_list
Expand All @@ -1046,6 +1095,18 @@
for p in o_block.phase_list
)

if mixed_block.include_inherent_reactions:
rhs += sum(
sum(
b.inherent_reaction_generation[t, o, p, j]
for j in mixed_block.component_list
if (p, j) in pc_set
)
for p in o_block.phase_list
)

return lhs == rhs

elif mb_type == MaterialBalanceType.elementTotal:
raise ConfigurationError(
"{} Separators do not support elemental "
Expand Down
Loading
Loading