#################################################################################
# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
# National Renewable Energy Laboratory, and National Energy Technology
# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
# of Energy). All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
# information, respectively. These files are also available online at the URL
# "https://github.com/watertap-org/watertap/"
#################################################################################
"""
ADM1 reaction package.
References:
[1] D.J. Batstone, J. Keller, I. Angelidaki, S.V. Kalyuzhnyi, S.G. Pavlostathis,
A. Rozzi, W.T.M. Sanders, H. Siegrist, V.A. Vavilin,
The IWA Anaerobic Digestion Model No 1 (ADM1), Water Science and Technology.
45 (2002) 65–73. https://doi.org/10.2166/wst.2002.0292.
[2] C. Rosén, U. Jeppsson, Aspects on ADM1 Implementation within the BSM2 Framework,
Department of Industrial Electrical Engineering and Automation, Lund University, Lund, Sweden. (2006) 1–35.
"""
# Import Pyomo libraries
import pyomo.environ as pyo
# Import IDAES cores
from idaes.core import (
declare_process_block_class,
MaterialFlowBasis,
ReactionParameterBlock,
ReactionBlockDataBase,
ReactionBlockBase,
)
from idaes.core.util.constants import Constants
from idaes.core.util.misc import add_object_reference
from idaes.core.util.exceptions import BurntToast
import idaes.logger as idaeslog
import idaes.core.util.scaling as iscale
from idaes.core.util.math import smooth_max
# Some more information about this module
__author__ = "Adam Atia, Alejandro Garciadiego, Xinhong Liu"
# Set up logger
_log = idaeslog.getLogger(__name__)
mw_n = 14 * pyo.units.kg / pyo.units.kmol
mw_c = 12 * pyo.units.kg / pyo.units.kmol
[docs]@declare_process_block_class("ADM1ReactionParameterBlock")
class ADM1ReactionParameterData(ReactionParameterBlock):
"""
Property Parameter Block Class
"""
[docs] def build(self):
"""
Callable method for Block construction.
"""
super().build()
self._reaction_block_class = ADM1ReactionBlock
# Reaction Index
# Reaction names based on standard numbering in ADM1
# R1: Disintegration
# R2: Hydrolysis of carbohydrates
# R3: Hydrolysis of proteins
# R4: Hydrolysis of lipids
# R5: Uptake of sugars
# R6: Uptake of amino acids
# R7: Uptake of long chain fatty acids (LCFAs)
# R8: Uptake of valerate
# R9: Uptake of butyrate
# R10: Uptake of propionate
# R11: Uptake of acetate
# R12: Uptake of hydrogen
# R13: Decay of X_su
# R14: Decay of X_aa
# R15: Decay of X_fa
# R16: Decay of X_c4
# R17: Decay of X_pro
# R18: Decay of X_ac
# R19: Decay of X_h2
self.rate_reaction_idx = pyo.Set(
initialize=[
"R1",
"R2",
"R3",
"R4",
"R5",
"R6",
"R7",
"R8",
"R9",
"R10",
"R11",
"R12",
"R13",
"R14",
"R15",
"R16",
"R17",
"R18",
"R19",
]
)
# Carbon content
Ci_dict = {
"S_su": 0.0313,
"S_aa": 0.03, # varies
"S_fa": 0.0217,
"S_va": 0.024,
"S_bu": 0.025,
"S_pro": 0.0268,
"S_ac": 0.0313,
"S_ch4": 0.0156,
"S_I": 0.03, # varies
"X_c": 0.02786, # varies
"X_ch": 0.0313,
"X_pr": 0.03, # varies
"X_li": 0.022,
"X_su": 0.0313,
"X_aa": 0.0313,
"X_fa": 0.0313,
"X_c4": 0.0313,
"X_pro": 0.0313,
"X_ac": 0.0313,
"X_h2": 0.0313,
"X_I": 0.03, # varies
}
self.Ci = pyo.Var(
Ci_dict.keys(),
initialize=Ci_dict,
units=pyo.units.kmol / pyo.units.kg,
domain=pyo.PositiveReals,
doc="Carbon content of component [kmole C/kg COD]",
)
# Stoichiometric Parameters (Table 6.1 in Batstone et al., 2002)
self.f_sI_xc = pyo.Var(
initialize=0.1,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Soluble inerts from composites",
)
self.f_xI_xc = pyo.Var(
initialize=0.20, # replacing 0.25 with 0.2 in accordance with Rosen & Jeppsson, 2006
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Particulate inerts from composites",
)
self.f_ch_xc = pyo.Var(
initialize=0.2,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Carbohydrates from composites",
)
self.f_pr_xc = pyo.Var(
initialize=0.2,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Proteins from composites",
)
self.f_li_xc = pyo.Var(
initialize=0.30, # replacing 0.25 with 0.3 in accordance with Rosen & Jeppsson, 2006
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Lipids from composites",
)
self.N_xc = pyo.Var(
initialize=0.0376
/ 14, # change from 0.002 to 0.0376/14 based on Rosen & Jeppsson, 2006
units=pyo.units.kmol * pyo.units.kg**-1,
domain=pyo.PositiveReals,
doc="Nitrogen content of composites [kmole N/kg COD]",
)
self.N_I = pyo.Var(
initialize=0.06
/ 14, # change from 0.002 to 0.06/14 based on Rosen & Jeppsson, 2006
units=pyo.units.kmol * pyo.units.kg**-1,
domain=pyo.PositiveReals,
doc="Nitrogen content of inerts [kmole N/kg COD]",
)
self.N_aa = pyo.Var(
initialize=0.007,
units=pyo.units.kmol * pyo.units.kg**-1,
domain=pyo.PositiveReals,
doc="Nitrogen in amino acids and proteins [kmole N/kg COD]",
)
self.N_bac = pyo.Var(
initialize=0.08 / 14,
units=pyo.units.kmol * pyo.units.kg**-1,
# units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Nitrogen content in bacteria [kmole N/kg COD]",
)
self.f_fa_li = pyo.Var(
initialize=0.95,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Fatty acids from lipids",
)
self.f_h2_su = pyo.Var(
initialize=0.19,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Hydrogen from sugars",
)
self.f_bu_su = pyo.Var(
initialize=0.13,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Butyrate from sugars",
)
self.f_pro_su = pyo.Var(
initialize=0.27,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Propionate from sugars",
)
self.f_ac_su = pyo.Var(
initialize=0.41,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Acetate from sugars",
)
self.f_h2_aa = pyo.Var(
initialize=0.06,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Hydrogen from amino acids",
)
self.f_va_aa = pyo.Var(
initialize=0.23,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Valerate from amino acids",
)
self.f_bu_aa = pyo.Var(
initialize=0.26,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Butyrate from amino acids",
)
self.f_pro_aa = pyo.Var(
initialize=0.05,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Propionate from amino acids",
)
self.f_ac_aa = pyo.Var(
initialize=0.40,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Acetate from amino acids",
)
self.Y_su = pyo.Var(
initialize=0.10,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on sugar substrate [kg COD X/ kg COD S]",
)
self.Y_aa = pyo.Var(
initialize=0.08,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on amino acid substrate [kg COD X/ kg COD S]",
)
self.Y_fa = pyo.Var(
initialize=0.06,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on fatty acid substrate [kg COD X/ kg COD S]",
)
self.Y_c4 = pyo.Var(
initialize=0.06,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on valerate and butyrate substrate [kg COD X/ kg COD S]",
)
self.Y_pro = pyo.Var(
initialize=0.04,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on propionate substrate [kg COD X/ kg COD S]",
)
self.Y_ac = pyo.Var(
initialize=0.05,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of biomass on acetate substrate [kg COD X/ kg COD S]",
)
self.Y_h2 = pyo.Var(
initialize=0.06,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Yield of hydrogen per biomass [kg COD S/ kg COD X]",
)
# Biochemical Parameters
self.k_dis = pyo.Var(
initialize=0.5,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order kinetic parameter for disintegration",
)
self.k_hyd_ch = pyo.Var(
initialize=10,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order kinetic parameter for hydrolysis of carbohydrates",
)
self.k_hyd_pr = pyo.Var(
initialize=10,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order kinetic parameter for hydrolysis of proteins",
)
self.k_hyd_li = pyo.Var(
initialize=10,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order kinetic parameter for hydrolysis of lipids",
)
self.K_S_IN = pyo.Var(
initialize=1e-4,
units=pyo.units.kmol * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Inhibition parameter for inorganic nitrogen",
)
self.k_m_su = pyo.Var(
initialize=30,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of sugars",
)
self.K_S_su = pyo.Var(
initialize=0.5,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of sugars",
)
self.pH_UL_aa = pyo.Var(
initialize=5.5,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Upper limit of pH for uptake rate of amino acids",
)
self.pH_LL_aa = pyo.Var(
initialize=4,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Lower limit of pH for uptake rate of amino acids",
)
self.k_m_aa = pyo.Var(
initialize=50,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of amino acids",
)
self.K_S_aa = pyo.Var(
initialize=0.3,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of amino acids",
)
self.k_m_fa = pyo.Var(
initialize=6,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of fatty acids",
)
self.K_S_fa = pyo.Var(
initialize=0.4,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of fatty acids",
)
self.K_I_h2_fa = pyo.Var(
initialize=5e-6,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Inhibition parameter for hydrogen during uptake of fatty acids",
)
self.k_m_c4 = pyo.Var(
initialize=20,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of valerate and butyrate",
)
self.K_S_c4 = pyo.Var(
initialize=0.2,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of valerate and butyrate",
)
self.K_I_h2_c4 = pyo.Var(
initialize=1e-5,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Inhibition parameter for hydrogen during uptake of valerate and butyrate",
)
self.k_m_pro = pyo.Var(
initialize=13,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of propionate",
)
self.K_S_pro = pyo.Var(
initialize=0.1,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of propionate",
)
self.K_I_h2_pro = pyo.Var(
initialize=3.5e-6,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Inhibition parameter for hydrogen during uptake of propionate",
)
self.k_m_ac = pyo.Var(
initialize=8,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of acetate",
)
self.K_S_ac = pyo.Var(
initialize=0.15,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of acetate",
)
self.K_I_nh3 = pyo.Var(
initialize=0.0018,
units=pyo.units.kmol * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Inhibition parameter for ammonia during uptake of acetate",
)
self.pH_UL_ac = pyo.Var(
initialize=7,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Upper limit of pH for uptake rate of acetate",
)
self.pH_LL_ac = pyo.Var(
initialize=6,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Lower limit of pH for uptake rate of acetate",
)
self.k_m_h2 = pyo.Var(
initialize=35,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="Monod maximum specific uptake rate of hydrogen",
)
self.K_S_h2 = pyo.Var(
initialize=7e-6,
units=pyo.units.kg * pyo.units.m**-3,
domain=pyo.PositiveReals,
doc="Half saturation value for uptake of hydrogen",
)
self.pH_UL_h2 = pyo.Var(
initialize=6,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Upper limit of pH for uptake rate of hydrogen",
)
self.pH_LL_h2 = pyo.Var(
initialize=5,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Lower limit of pH for uptake rate of hydrogen",
)
self.k_dec_X_su = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_su",
)
self.k_dec_X_aa = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_aa",
)
self.k_dec_X_fa = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_fa",
)
self.k_dec_X_c4 = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_c4",
)
self.k_dec_X_pro = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_pro",
)
self.k_dec_X_ac = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_ac",
)
self.k_dec_X_h2 = pyo.Var(
initialize=0.02,
units=pyo.units.day**-1,
domain=pyo.PositiveReals,
doc="First-order decay rate for X_h2",
)
self.K_a_va = pyo.Var(
initialize=10 ** (-4.86),
units=pyo.units.kmol / pyo.units.m**3,
domain=pyo.PositiveReals,
doc="Valerate acid-base equilibrium constant",
)
self.K_a_bu = pyo.Var(
initialize=10 ** (-4.82),
units=pyo.units.kmol / pyo.units.m**3,
domain=pyo.PositiveReals,
doc="Butyrate acid-base equilibrium constant",
)
self.K_a_pro = pyo.Var(
initialize=10 ** (-4.88),
units=pyo.units.kmol / pyo.units.m**3,
domain=pyo.PositiveReals,
doc="Propionate acid-base equilibrium constant",
)
self.K_a_ac = pyo.Var(
initialize=10 ** (-4.76),
units=pyo.units.kmol / pyo.units.m**3,
domain=pyo.PositiveReals,
doc="Acetate acid-base equilibrium constant",
)
self.temperature_ref = pyo.Param(
within=pyo.PositiveReals,
mutable=True,
default=298.15,
doc="Reference temperature",
units=pyo.units.K,
)
# Reaction Stoichiometry
# This is the stoichiometric part of the Peterson matrix in dict form.
# See Table 3.1 in Batstone et al., 2002.
# Exclude non-zero stoichiometric coefficients for S_IC initially since they depend on other stoichiometric coefficients.
self.rate_reaction_stoichiometry = {
# R1: Disintegration
("R1", "Liq", "H2O"): 0,
("R1", "Liq", "S_su"): 0,
("R1", "Liq", "S_aa"): 0,
("R1", "Liq", "S_fa"): 0,
("R1", "Liq", "S_va"): 0,
("R1", "Liq", "S_bu"): 0,
("R1", "Liq", "S_pro"): 0,
("R1", "Liq", "S_ac"): 0,
("R1", "Liq", "S_h2"): 0,
("R1", "Liq", "S_ch4"): 0,
("R1", "Liq", "S_IC"): (
self.Ci["X_c"]
- self.f_sI_xc * self.Ci["S_I"]
- self.f_ch_xc * self.Ci["X_ch"]
- self.f_pr_xc * self.Ci["X_pr"]
- self.f_li_xc * self.Ci["X_li"]
- self.f_xI_xc * self.Ci["X_I"]
)
* mw_c,
("R1", "Liq", "S_IN"): (
self.N_xc
- self.f_xI_xc * self.N_I
- self.f_sI_xc * self.N_I
- self.f_pr_xc * self.N_aa
)
* mw_n,
("R1", "Liq", "S_I"): self.f_sI_xc,
("R1", "Liq", "X_c"): -1,
("R1", "Liq", "X_ch"): self.f_ch_xc,
("R1", "Liq", "X_pr"): self.f_pr_xc,
("R1", "Liq", "X_li"): self.f_li_xc,
("R1", "Liq", "X_su"): 0,
("R1", "Liq", "X_aa"): 0,
("R1", "Liq", "X_fa"): 0,
("R1", "Liq", "X_c4"): 0,
("R1", "Liq", "X_pro"): 0,
("R1", "Liq", "X_ac"): 0,
("R1", "Liq", "X_h2"): 0,
("R1", "Liq", "X_I"): self.f_xI_xc,
# R2: Hydrolysis of carbohydrates
("R2", "Liq", "H2O"): 0,
("R2", "Liq", "S_su"): 1,
("R2", "Liq", "S_aa"): 0,
("R2", "Liq", "S_fa"): 0,
("R2", "Liq", "S_va"): 0,
("R2", "Liq", "S_bu"): 0,
("R2", "Liq", "S_pro"): 0,
("R2", "Liq", "S_ac"): 0,
("R2", "Liq", "S_h2"): 0,
("R2", "Liq", "S_ch4"): 0,
("R2", "Liq", "S_IC"): 0,
("R2", "Liq", "S_IN"): 0,
("R2", "Liq", "S_I"): 0,
("R2", "Liq", "X_c"): 0,
("R2", "Liq", "X_ch"): -1,
("R2", "Liq", "X_pr"): 0,
("R2", "Liq", "X_li"): 0,
("R2", "Liq", "X_su"): 0,
("R2", "Liq", "X_aa"): 0,
("R2", "Liq", "X_fa"): 0,
("R2", "Liq", "X_c4"): 0,
("R2", "Liq", "X_pro"): 0,
("R2", "Liq", "X_ac"): 0,
("R2", "Liq", "X_h2"): 0,
("R2", "Liq", "X_I"): 0,
# R3: Hydrolysis of proteins
("R3", "Liq", "H2O"): 0,
("R3", "Liq", "S_su"): 0,
("R3", "Liq", "S_aa"): 1,
("R3", "Liq", "S_fa"): 0,
("R3", "Liq", "S_va"): 0,
("R3", "Liq", "S_bu"): 0,
("R3", "Liq", "S_pro"): 0,
("R3", "Liq", "S_ac"): 0,
("R3", "Liq", "S_h2"): 0,
("R3", "Liq", "S_ch4"): 0,
("R3", "Liq", "S_IC"): 0,
("R3", "Liq", "S_IN"): 0,
("R3", "Liq", "S_I"): 0,
("R3", "Liq", "X_c"): 0,
("R3", "Liq", "X_ch"): 0,
("R3", "Liq", "X_pr"): -1,
("R3", "Liq", "X_li"): 0,
("R3", "Liq", "X_su"): 0,
("R3", "Liq", "X_aa"): 0,
("R3", "Liq", "X_fa"): 0,
("R3", "Liq", "X_c4"): 0,
("R3", "Liq", "X_pro"): 0,
("R3", "Liq", "X_ac"): 0,
("R3", "Liq", "X_h2"): 0,
("R3", "Liq", "X_I"): 0,
# R4: Hydrolysis of lipids
("R4", "Liq", "H2O"): 0,
("R4", "Liq", "S_su"): 1 - self.f_fa_li,
("R4", "Liq", "S_aa"): 0,
("R4", "Liq", "S_fa"): self.f_fa_li,
("R4", "Liq", "S_va"): 0,
("R4", "Liq", "S_bu"): 0,
("R4", "Liq", "S_pro"): 0,
("R4", "Liq", "S_ac"): 0,
("R4", "Liq", "S_h2"): 0,
("R4", "Liq", "S_ch4"): 0,
("R4", "Liq", "S_IC"): (
-(1 - self.f_fa_li) * self.Ci["S_su"]
- self.f_fa_li * self.Ci["S_fa"]
+ self.Ci["X_li"]
)
* mw_c,
("R4", "Liq", "S_IN"): 0,
("R4", "Liq", "S_I"): 0,
("R4", "Liq", "X_c"): 0,
("R4", "Liq", "X_ch"): 0,
("R4", "Liq", "X_pr"): 0,
("R4", "Liq", "X_li"): -1,
("R4", "Liq", "X_su"): 0,
("R4", "Liq", "X_aa"): 0,
("R4", "Liq", "X_fa"): 0,
("R4", "Liq", "X_c4"): 0,
("R4", "Liq", "X_pro"): 0,
("R4", "Liq", "X_ac"): 0,
("R4", "Liq", "X_h2"): 0,
("R4", "Liq", "X_I"): 0,
# R5: Uptake of sugars
("R5", "Liq", "H2O"): 0,
("R5", "Liq", "S_su"): -1,
("R5", "Liq", "S_aa"): 0,
("R5", "Liq", "S_fa"): 0,
("R5", "Liq", "S_va"): 0,
("R5", "Liq", "S_bu"): (1 - self.Y_su) * self.f_bu_su,
("R5", "Liq", "S_pro"): (1 - self.Y_su) * self.f_pro_su,
("R5", "Liq", "S_ac"): (1 - self.Y_su) * self.f_ac_su,
("R5", "Liq", "S_h2"): (1 - self.Y_su) * self.f_h2_su,
("R5", "Liq", "S_ch4"): 0,
("R5", "Liq", "S_IC"): (
self.Ci["S_su"]
- (1 - self.Y_su)
* (
self.f_bu_su * self.Ci["S_bu"]
+ self.f_pro_su * self.Ci["S_pro"]
+ self.f_ac_su * self.Ci["S_ac"]
)
- self.Y_su * self.Ci["X_su"]
)
* mw_c,
("R5", "Liq", "S_IN"): (-self.Y_su * self.N_bac) * mw_n,
("R5", "Liq", "S_I"): 0,
("R5", "Liq", "X_c"): 0,
("R5", "Liq", "X_ch"): 0,
("R5", "Liq", "X_pr"): 0,
("R5", "Liq", "X_li"): 0,
("R5", "Liq", "X_su"): self.Y_su,
("R5", "Liq", "X_aa"): 0,
("R5", "Liq", "X_fa"): 0,
("R5", "Liq", "X_c4"): 0,
("R5", "Liq", "X_pro"): 0,
("R5", "Liq", "X_ac"): 0,
("R5", "Liq", "X_h2"): 0,
("R5", "Liq", "X_I"): 0,
# R6: Uptake of amino acids
("R6", "Liq", "H2O"): 0,
("R6", "Liq", "S_su"): 0,
("R6", "Liq", "S_aa"): -1,
("R6", "Liq", "S_fa"): 0,
("R6", "Liq", "S_va"): (1 - self.Y_aa) * self.f_va_aa,
("R6", "Liq", "S_bu"): (1 - self.Y_aa) * self.f_bu_aa,
("R6", "Liq", "S_pro"): (1 - self.Y_aa) * self.f_pro_aa,
("R6", "Liq", "S_ac"): (1 - self.Y_aa) * self.f_ac_aa,
("R6", "Liq", "S_h2"): (1 - self.Y_aa) * self.f_h2_aa,
("R6", "Liq", "S_ch4"): 0,
("R6", "Liq", "S_IC"): (
self.Ci["S_aa"]
- (1 - self.Y_aa)
* (
self.f_va_aa * self.Ci["S_va"]
+ self.f_bu_aa * self.Ci["S_bu"]
+ self.f_pro_aa * self.Ci["S_pro"]
+ self.f_ac_aa * self.Ci["S_ac"]
)
- self.Y_aa * self.Ci["X_aa"]
)
* mw_c,
("R6", "Liq", "S_IN"): (self.N_aa - self.Y_aa * self.N_bac) * mw_n,
("R6", "Liq", "S_I"): 0,
("R6", "Liq", "X_c"): 0,
("R6", "Liq", "X_ch"): 0,
("R6", "Liq", "X_pr"): 0,
("R6", "Liq", "X_li"): 0,
("R6", "Liq", "X_su"): 0,
("R6", "Liq", "X_aa"): self.Y_aa,
("R6", "Liq", "X_fa"): 0,
("R6", "Liq", "X_c4"): 0,
("R6", "Liq", "X_pro"): 0,
("R6", "Liq", "X_ac"): 0,
("R6", "Liq", "X_h2"): 0,
("R6", "Liq", "X_I"): 0,
# R7: Uptake of long chain fatty acids (LCFAs)
("R7", "Liq", "H2O"): 0,
("R7", "Liq", "S_su"): 0,
("R7", "Liq", "S_aa"): 0,
("R7", "Liq", "S_fa"): -1,
("R7", "Liq", "S_va"): 0,
("R7", "Liq", "S_bu"): 0,
("R7", "Liq", "S_pro"): 0,
("R7", "Liq", "S_ac"): (1 - self.Y_fa) * 0.7,
("R7", "Liq", "S_h2"): (1 - self.Y_fa) * 0.3,
("R7", "Liq", "S_ch4"): 0,
("R7", "Liq", "S_IC"): (
self.Ci["S_fa"]
- (1 - self.Y_fa) * 0.7 * self.Ci["S_ac"]
- self.Y_fa * self.Ci["X_fa"]
)
* mw_c,
("R7", "Liq", "S_IN"): (-self.Y_fa * self.N_bac) * mw_n,
("R7", "Liq", "S_I"): 0,
("R7", "Liq", "X_c"): 0,
("R7", "Liq", "X_ch"): 0,
("R7", "Liq", "X_pr"): 0,
("R7", "Liq", "X_li"): 0,
("R7", "Liq", "X_su"): 0,
("R7", "Liq", "X_aa"): 0,
("R7", "Liq", "X_fa"): self.Y_fa,
("R7", "Liq", "X_c4"): 0,
("R7", "Liq", "X_pro"): 0,
("R7", "Liq", "X_ac"): 0,
("R7", "Liq", "X_h2"): 0,
("R7", "Liq", "X_I"): 0,
# R8: Uptake of valerate
("R8", "Liq", "H2O"): 0,
("R8", "Liq", "S_su"): 0,
("R8", "Liq", "S_aa"): 0,
("R8", "Liq", "S_fa"): 0,
("R8", "Liq", "S_va"): -1,
("R8", "Liq", "S_bu"): 0,
("R8", "Liq", "S_pro"): (1 - self.Y_c4) * 0.54,
("R8", "Liq", "S_ac"): (1 - self.Y_c4) * 0.31,
("R8", "Liq", "S_h2"): (1 - self.Y_c4) * 0.15,
("R8", "Liq", "S_ch4"): 0,
("R8", "Liq", "S_IC"): (
self.Ci["S_va"]
- (1 - self.Y_c4) * 0.54 * self.Ci["S_pro"]
- (1 - self.Y_c4) * 0.31 * self.Ci["S_ac"]
- self.Y_c4 * self.Ci["X_c4"]
)
* mw_c,
("R8", "Liq", "S_IN"): (-self.Y_c4 * self.N_bac) * mw_n,
("R8", "Liq", "S_I"): 0,
("R8", "Liq", "X_c"): 0,
("R8", "Liq", "X_ch"): 0,
("R8", "Liq", "X_pr"): 0,
("R8", "Liq", "X_li"): 0,
("R8", "Liq", "X_su"): 0,
("R8", "Liq", "X_aa"): 0,
("R8", "Liq", "X_fa"): 0,
("R8", "Liq", "X_c4"): self.Y_c4,
("R8", "Liq", "X_pro"): 0,
("R8", "Liq", "X_ac"): 0,
("R8", "Liq", "X_h2"): 0,
("R8", "Liq", "X_I"): 0,
# R9: Uptake of butyrate
("R9", "Liq", "H2O"): 0,
("R9", "Liq", "S_su"): 0,
("R9", "Liq", "S_aa"): 0,
("R9", "Liq", "S_fa"): 0,
("R9", "Liq", "S_va"): 0,
("R9", "Liq", "S_bu"): -1,
("R9", "Liq", "S_pro"): 0,
("R9", "Liq", "S_ac"): (1 - self.Y_c4) * 0.8,
("R9", "Liq", "S_h2"): (1 - self.Y_c4) * 0.2,
("R9", "Liq", "S_ch4"): 0,
("R9", "Liq", "S_IC"): (
self.Ci["S_bu"]
- (1 - self.Y_c4) * 0.8 * self.Ci["S_ac"]
- self.Y_c4 * self.Ci["X_c4"]
)
* mw_c,
("R9", "Liq", "S_IN"): (-self.Y_c4 * self.N_bac) * mw_n,
("R9", "Liq", "S_I"): 0,
("R9", "Liq", "X_c"): 0,
("R9", "Liq", "X_ch"): 0,
("R9", "Liq", "X_pr"): 0,
("R9", "Liq", "X_li"): 0,
("R9", "Liq", "X_su"): 0,
("R9", "Liq", "X_aa"): 0,
("R9", "Liq", "X_fa"): 0,
("R9", "Liq", "X_c4"): self.Y_c4,
("R9", "Liq", "X_pro"): 0,
("R9", "Liq", "X_ac"): 0,
("R9", "Liq", "X_h2"): 0,
("R9", "Liq", "X_I"): 0,
# R10: Uptake of propionate
("R10", "Liq", "H2O"): 0,
("R10", "Liq", "S_su"): 0,
("R10", "Liq", "S_aa"): 0,
("R10", "Liq", "S_fa"): 0,
("R10", "Liq", "S_va"): 0,
("R10", "Liq", "S_bu"): 0,
("R10", "Liq", "S_pro"): -1,
("R10", "Liq", "S_ac"): (1 - self.Y_pro) * 0.57,
("R10", "Liq", "S_h2"): (1 - self.Y_pro) * 0.43,
("R10", "Liq", "S_ch4"): 0,
("R10", "Liq", "S_IC"): (
self.Ci["S_pro"]
- (1 - self.Y_pro) * 0.57 * self.Ci["S_ac"]
- self.Y_pro * self.Ci["X_pro"]
)
* mw_c,
("R10", "Liq", "S_IN"): (-self.Y_pro * self.N_bac) * mw_n,
("R10", "Liq", "S_I"): 0,
("R10", "Liq", "X_c"): 0,
("R10", "Liq", "X_ch"): 0,
("R10", "Liq", "X_pr"): 0,
("R10", "Liq", "X_li"): 0,
("R10", "Liq", "X_su"): 0,
("R10", "Liq", "X_aa"): 0,
("R10", "Liq", "X_fa"): 0,
("R10", "Liq", "X_c4"): 0,
("R10", "Liq", "X_pro"): self.Y_pro,
("R10", "Liq", "X_ac"): 0,
("R10", "Liq", "X_h2"): 0,
("R10", "Liq", "X_I"): 0,
# R11: Uptake of acetate
("R11", "Liq", "H2O"): 0,
("R11", "Liq", "S_su"): 0,
("R11", "Liq", "S_aa"): 0,
("R11", "Liq", "S_fa"): 0,
("R11", "Liq", "S_va"): 0,
("R11", "Liq", "S_bu"): 0,
("R11", "Liq", "S_pro"): 0,
("R11", "Liq", "S_ac"): -1,
("R11", "Liq", "S_h2"): 0,
("R11", "Liq", "S_ch4"): 1 - self.Y_ac,
("R11", "Liq", "S_IC"): (
self.Ci["S_ac"]
- (1 - self.Y_ac) * self.Ci["S_ch4"]
- self.Y_ac * self.Ci["X_ac"]
)
* mw_c,
("R11", "Liq", "S_IN"): (-self.Y_ac * self.N_bac) * mw_n,
("R11", "Liq", "S_I"): 0,
("R11", "Liq", "X_c"): 0,
("R11", "Liq", "X_ch"): 0,
("R11", "Liq", "X_pr"): 0,
("R11", "Liq", "X_li"): 0,
("R11", "Liq", "X_su"): 0,
("R11", "Liq", "X_aa"): 0,
("R11", "Liq", "X_fa"): 0,
("R11", "Liq", "X_c4"): 0,
("R11", "Liq", "X_pro"): 0,
("R11", "Liq", "X_ac"): self.Y_ac,
("R11", "Liq", "X_h2"): 0,
("R11", "Liq", "X_I"): 0,
# R12: Uptake of hydrogen
("R12", "Liq", "H2O"): 0,
("R12", "Liq", "S_su"): 0,
("R12", "Liq", "S_aa"): 0,
("R12", "Liq", "S_fa"): 0,
("R12", "Liq", "S_va"): 0,
("R12", "Liq", "S_bu"): 0,
("R12", "Liq", "S_pro"): 0,
("R12", "Liq", "S_ac"): 0,
("R12", "Liq", "S_h2"): -1,
("R12", "Liq", "S_ch4"): (1 - self.Y_h2),
("R12", "Liq", "S_IC"): (
-(1 - self.Y_h2) * self.Ci["S_ch4"] - self.Y_h2 * self.Ci["X_h2"]
)
* mw_c,
("R12", "Liq", "S_IN"): (-self.Y_h2 * self.N_bac) * mw_n,
("R12", "Liq", "S_I"): 0,
("R12", "Liq", "X_c"): 0,
("R12", "Liq", "X_ch"): 0,
("R12", "Liq", "X_pr"): 0,
("R12", "Liq", "X_li"): 0,
("R12", "Liq", "X_su"): 0,
("R12", "Liq", "X_aa"): 0,
("R12", "Liq", "X_fa"): 0,
("R12", "Liq", "X_c4"): 0,
("R12", "Liq", "X_pro"): 0,
("R12", "Liq", "X_ac"): 0,
("R12", "Liq", "X_h2"): self.Y_h2,
("R12", "Liq", "X_I"): 0,
# R13: Decay of X_su
("R13", "Liq", "H2O"): 0,
("R13", "Liq", "S_su"): 0,
("R13", "Liq", "S_aa"): 0,
("R13", "Liq", "S_fa"): 0,
("R13", "Liq", "S_va"): 0,
("R13", "Liq", "S_bu"): 0,
("R13", "Liq", "S_pro"): 0,
("R13", "Liq", "S_ac"): 0,
("R13", "Liq", "S_h2"): 0,
("R13", "Liq", "S_ch4"): 0,
("R13", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R13", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R13", "Liq", "S_I"): 0,
("R13", "Liq", "X_c"): 1,
("R13", "Liq", "X_ch"): 0,
("R13", "Liq", "X_pr"): 0,
("R13", "Liq", "X_li"): 0,
("R13", "Liq", "X_su"): -1,
("R13", "Liq", "X_aa"): 0,
("R13", "Liq", "X_fa"): 0,
("R13", "Liq", "X_c4"): 0,
("R13", "Liq", "X_pro"): 0,
("R13", "Liq", "X_ac"): 0,
("R13", "Liq", "X_h2"): 0,
("R13", "Liq", "X_I"): 0,
# R14: Decay of X_aa
("R14", "Liq", "H2O"): 0,
("R14", "Liq", "S_su"): 0,
("R14", "Liq", "S_aa"): 0,
("R14", "Liq", "S_fa"): 0,
("R14", "Liq", "S_va"): 0,
("R14", "Liq", "S_bu"): 0,
("R14", "Liq", "S_pro"): 0,
("R14", "Liq", "S_ac"): 0,
("R14", "Liq", "S_h2"): 0,
("R14", "Liq", "S_ch4"): 0,
("R14", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R14", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R14", "Liq", "S_I"): 0,
("R14", "Liq", "X_c"): 1,
("R14", "Liq", "X_ch"): 0,
("R14", "Liq", "X_pr"): 0,
("R14", "Liq", "X_li"): 0,
("R14", "Liq", "X_su"): 0,
("R14", "Liq", "X_aa"): -1,
("R14", "Liq", "X_fa"): 0,
("R14", "Liq", "X_c4"): 0,
("R14", "Liq", "X_pro"): 0,
("R14", "Liq", "X_ac"): 0,
("R14", "Liq", "X_h2"): 0,
("R14", "Liq", "X_I"): 0,
# R15: Decay of X_fa
("R15", "Liq", "H2O"): 0,
("R15", "Liq", "S_su"): 0,
("R15", "Liq", "S_aa"): 0,
("R15", "Liq", "S_fa"): 0,
("R15", "Liq", "S_va"): 0,
("R15", "Liq", "S_bu"): 0,
("R15", "Liq", "S_pro"): 0,
("R15", "Liq", "S_ac"): 0,
("R15", "Liq", "S_h2"): 0,
("R15", "Liq", "S_ch4"): 0,
("R15", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R15", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R15", "Liq", "S_I"): 0,
("R15", "Liq", "X_c"): 1,
("R15", "Liq", "X_ch"): 0,
("R15", "Liq", "X_pr"): 0,
("R15", "Liq", "X_li"): 0,
("R15", "Liq", "X_su"): 0,
("R15", "Liq", "X_aa"): 0,
("R15", "Liq", "X_fa"): -1,
("R15", "Liq", "X_c4"): 0,
("R15", "Liq", "X_pro"): 0,
("R15", "Liq", "X_ac"): 0,
("R15", "Liq", "X_h2"): 0,
("R15", "Liq", "X_I"): 0,
# R16: Decay of X_c4
("R16", "Liq", "H2O"): 0,
("R16", "Liq", "S_su"): 0,
("R16", "Liq", "S_aa"): 0,
("R16", "Liq", "S_fa"): 0,
("R16", "Liq", "S_va"): 0,
("R16", "Liq", "S_bu"): 0,
("R16", "Liq", "S_pro"): 0,
("R16", "Liq", "S_ac"): 0,
("R16", "Liq", "S_h2"): 0,
("R16", "Liq", "S_ch4"): 0,
("R16", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R16", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R16", "Liq", "S_I"): 0,
("R16", "Liq", "X_c"): 1,
("R16", "Liq", "X_ch"): 0,
("R16", "Liq", "X_pr"): 0,
("R16", "Liq", "X_li"): 0,
("R16", "Liq", "X_su"): 0,
("R16", "Liq", "X_aa"): 0,
("R16", "Liq", "X_fa"): 0,
("R16", "Liq", "X_c4"): -1,
("R16", "Liq", "X_pro"): 0,
("R16", "Liq", "X_ac"): 0,
("R16", "Liq", "X_h2"): 0,
("R16", "Liq", "X_I"): 0,
# R17: Decay of X_pro
("R17", "Liq", "H2O"): 0,
("R17", "Liq", "S_su"): 0,
("R17", "Liq", "S_aa"): 0,
("R17", "Liq", "S_fa"): 0,
("R17", "Liq", "S_va"): 0,
("R17", "Liq", "S_bu"): 0,
("R17", "Liq", "S_pro"): 0,
("R17", "Liq", "S_ac"): 0,
("R17", "Liq", "S_h2"): 0,
("R17", "Liq", "S_ch4"): 0,
("R17", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R17", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R17", "Liq", "S_I"): 0,
("R17", "Liq", "X_c"): 1,
("R17", "Liq", "X_ch"): 0,
("R17", "Liq", "X_pr"): 0,
("R17", "Liq", "X_li"): 0,
("R17", "Liq", "X_su"): 0,
("R17", "Liq", "X_aa"): 0,
("R17", "Liq", "X_fa"): 0,
("R17", "Liq", "X_c4"): 0,
("R17", "Liq", "X_pro"): -1,
("R17", "Liq", "X_ac"): 0,
("R17", "Liq", "X_h2"): 0,
("R17", "Liq", "X_I"): 0,
# R18: Decay of X_ac
("R18", "Liq", "H2O"): 0,
("R18", "Liq", "S_su"): 0,
("R18", "Liq", "S_aa"): 0,
("R18", "Liq", "S_fa"): 0,
("R18", "Liq", "S_va"): 0,
("R18", "Liq", "S_bu"): 0,
("R18", "Liq", "S_pro"): 0,
("R18", "Liq", "S_ac"): 0,
("R18", "Liq", "S_h2"): 0,
("R18", "Liq", "S_ch4"): 0,
("R18", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R18", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R18", "Liq", "S_I"): 0,
("R18", "Liq", "X_c"): 1,
("R18", "Liq", "X_ch"): 0,
("R18", "Liq", "X_pr"): 0,
("R18", "Liq", "X_li"): 0,
("R18", "Liq", "X_su"): 0,
("R18", "Liq", "X_aa"): 0,
("R18", "Liq", "X_fa"): 0,
("R18", "Liq", "X_c4"): 0,
("R18", "Liq", "X_pro"): 0,
("R18", "Liq", "X_ac"): -1,
("R18", "Liq", "X_h2"): 0,
("R18", "Liq", "X_I"): 0,
# R19: Decay of X_h2
("R19", "Liq", "H2O"): 0,
("R19", "Liq", "S_su"): 0,
("R19", "Liq", "S_aa"): 0,
("R19", "Liq", "S_fa"): 0,
("R19", "Liq", "S_va"): 0,
("R19", "Liq", "S_bu"): 0,
("R19", "Liq", "S_pro"): 0,
("R19", "Liq", "S_ac"): 0,
("R19", "Liq", "S_h2"): 0,
("R19", "Liq", "S_ch4"): 0,
("R19", "Liq", "S_IC"): (-self.Ci["X_c"] + self.Ci["X_ac"]) * mw_c,
("R19", "Liq", "S_IN"): (self.N_bac - self.N_xc) * mw_n,
("R19", "Liq", "S_I"): 0,
("R19", "Liq", "X_c"): 1,
("R19", "Liq", "X_ch"): 0,
("R19", "Liq", "X_pr"): 0,
("R19", "Liq", "X_li"): 0,
("R19", "Liq", "X_su"): 0,
("R19", "Liq", "X_aa"): 0,
("R19", "Liq", "X_fa"): 0,
("R19", "Liq", "X_c4"): 0,
("R19", "Liq", "X_pro"): 0,
("R19", "Liq", "X_ac"): 0,
("R19", "Liq", "X_h2"): -1,
("R19", "Liq", "X_I"): 0,
}
for R in self.rate_reaction_idx:
self.rate_reaction_stoichiometry[R, "Liq", "S_cat"] = 0
self.rate_reaction_stoichiometry[R, "Liq", "S_an"] = 0
# Fix all the variables we just created
for v in self.component_objects(pyo.Var, descend_into=False):
v.fix()
[docs]class _ADM1ReactionBlock(ReactionBlockBase):
"""
This Class contains methods which should be applied to Reaction Blocks as a
whole, rather than individual elements of indexed Reaction Blocks.
"""
[docs] def initialize(self, outlvl=idaeslog.NOTSET, **kwargs):
"""
Initialization routine for reaction package.
Keyword Arguments:
outlvl : sets output level of initialization routine
Returns:
None
"""
init_log = idaeslog.getInitLogger(self.name, outlvl, tag="properties")
init_log.info("Initialization Complete.")
[docs]@declare_process_block_class("ADM1ReactionBlock", block_class=_ADM1ReactionBlock)
class ADM1ReactionBlockData(ReactionBlockDataBase):
"""
ReactionBlock for ADM1.
"""
[docs] def build(self):
"""
Callable method for Block construction
"""
super().build()
# Create references to state vars
# Concentration
add_object_reference(self, "conc_mass_comp_ref", self.state_ref.conc_mass_comp)
add_object_reference(self, "temperature", self.state_ref.temperature)
# Initial values of rates of reaction [2]
self.rates = {
"R1": 1.786e-06,
"R2": 3.235e-06,
"R3": 1.187e-05,
"R4": 3.412e-06,
"R5": 3.404e-06,
"R6": 1.187e-05,
"R7": 3.185e-06,
"R8": 2.505e-06,
"R9": 3.230e-06,
"R10": 2.636e-06,
"R11": 1.220e-05,
"R12": 4.184e-06,
"R13": 9.736e-08,
"R14": 2.730e-07,
"R15": 5.627e-08,
"R16": 9.999e-08,
"R17": 3.179e-08,
"R18": 1.761e-07,
"R19": 7.339e-08,
}
# Rate of reaction method
def _rxn_rate(self):
self.reaction_rate = pyo.Var(
self.params.rate_reaction_idx,
initialize=self.rates,
domain=pyo.NonNegativeReals,
doc="Rate of reaction",
units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
self.I = pyo.Var(
self.params.rate_reaction_idx,
initialize=1,
bounds=(1e-8, 10),
doc="Process inhibition term",
units=pyo.units.dimensionless,
)
self.pKW = pyo.Var(
initialize=14,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Water dissociation constant",
)
self.pK_a_co2 = pyo.Var(
initialize=6.35,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Carbon dioxide acid-base equilibrium constant",
)
self.pK_a_IN = pyo.Var(
initialize=9.25,
units=pyo.units.dimensionless,
domain=pyo.PositiveReals,
doc="Inorganic nitrogen acid-base equilibrium constant",
)
self.conc_mass_va = pyo.Var(
initialize=0.011,
domain=pyo.NonNegativeReals,
doc="mass concentration of va-",
units=pyo.units.kg / pyo.units.m**3,
)
self.conc_mass_bu = pyo.Var(
initialize=0.013,
domain=pyo.NonNegativeReals,
doc="mass concentration of bu-",
units=pyo.units.kg / pyo.units.m**3,
)
self.conc_mass_pro = pyo.Var(
initialize=0.016,
domain=pyo.NonNegativeReals,
doc="mass concentration of pro-",
units=pyo.units.kg / pyo.units.m**3,
)
self.conc_mass_ac = pyo.Var(
initialize=0.2,
domain=pyo.NonNegativeReals,
doc="mass concentration of ac-",
units=pyo.units.kg / pyo.units.m**3,
)
self.conc_mol_hco3 = pyo.Var(
initialize=0.14,
domain=pyo.NonNegativeReals,
doc="molar concentration of hco3",
units=pyo.units.kmol / pyo.units.m**3,
)
self.conc_mol_nh3 = pyo.Var(
initialize=0.0041,
domain=pyo.NonNegativeReals,
doc="molar concentration of nh3",
units=pyo.units.kmol / pyo.units.m**3,
)
self.conc_mol_co2 = pyo.Var(
initialize=0.0099,
domain=pyo.NonNegativeReals,
doc="molar concentration of co2",
units=pyo.units.kmol / pyo.units.m**3,
)
self.conc_mol_nh4 = pyo.Var(
initialize=0.1261,
domain=pyo.NonNegativeReals,
doc="molar concentration of nh4",
units=pyo.units.kmol / pyo.units.m**3,
)
self.S_H = pyo.Var(
initialize=1e-7,
domain=pyo.NonNegativeReals,
doc="molar concentration of H",
units=pyo.units.kmol / pyo.units.m**3,
)
self.pH = pyo.Var(
initialize=7,
bounds=(0, 14),
domain=pyo.PositiveReals,
doc="pH of solution",
units=pyo.units.dimensionless,
)
# Equation from [2]
def Dissociation_rule(self, t):
return pyo.log(10**-self.pKW) == (
pyo.log(1e-14)
+ 55900
/ pyo.units.mole
* pyo.units.joule
/ (Constants.gas_constant)
* ((1 / self.params.temperature_ref) - (1 / self.temperature))
)
self.Dissociation = pyo.Constraint(
rule=Dissociation_rule,
doc="Water dissociation constant constraint",
)
# Equation from [2]
def CO2_acid_base_equilibrium_rule(self, t):
return pyo.log(10**-self.pK_a_co2) == (
pyo.log(10**-6.35)
+ 7646
/ pyo.units.mole
* pyo.units.joule
/ (Constants.gas_constant)
* ((1 / self.params.temperature_ref) - (1 / self.temperature))
)
self.CO2_acid_base_equilibrium = pyo.Constraint(
rule=CO2_acid_base_equilibrium_rule,
doc="Carbon dioxide acid-base equilibrium constraint",
)
# Equation from [2]
def IN_acid_base_equilibrium_rule(self, t):
return pyo.log(10**-self.pK_a_IN) == (
pyo.log(10**-9.25)
+ 51965
/ pyo.units.mole
* pyo.units.joule
/ (Constants.gas_constant)
* ((1 / self.params.temperature_ref) - (1 / self.temperature))
)
self.IN_acid_base_equilibrium = pyo.Constraint(
rule=IN_acid_base_equilibrium_rule,
doc="Nitrogen acid-base equilibrium constraint",
)
def rule_pH(self):
return self.pH == -pyo.log10(self.S_H / (pyo.units.kmole / pyo.units.m**3))
self.pH_calc = pyo.Constraint(rule=rule_pH, doc="pH of solution")
def concentration_of_va_rule(self):
return (
self.conc_mass_va * (1 + self.S_H / self.params.K_a_va)
== self.conc_mass_comp_ref["S_va"]
)
self.concentration_of_va = pyo.Constraint(
rule=concentration_of_va_rule,
doc="constraint concentration of va-",
)
def concentration_of_bu_rule(self):
return (
self.conc_mass_bu * (1 + self.S_H / self.params.K_a_bu)
== self.conc_mass_comp_ref["S_bu"]
)
self.concentration_of_bu = pyo.Constraint(
rule=concentration_of_bu_rule,
doc="constraint concentration of bu-",
)
def concentration_of_pro_rule(self):
return (
self.conc_mass_pro * (1 + self.S_H / self.params.K_a_pro)
== self.conc_mass_comp_ref["S_pro"]
)
self.concentration_of_pro = pyo.Constraint(
rule=concentration_of_pro_rule,
doc="constraint concentration of pro-",
)
def concentration_of_ac_rule(self):
return (
self.conc_mass_ac * (1 + self.S_H / self.params.K_a_ac)
== self.conc_mass_comp_ref["S_ac"]
)
self.concentration_of_ac = pyo.Constraint(
rule=concentration_of_ac_rule,
doc="constraint concentration of ac-",
)
def concentration_of_hco3_rule(self):
return (
self.pK_a_co2
== pyo.log10(self.conc_mol_co2 / (pyo.units.kmole / pyo.units.m**3))
- pyo.log10(self.conc_mol_hco3 / (pyo.units.kmole / pyo.units.m**3))
+ self.pH
)
self.concentration_of_hco3 = pyo.Constraint(
rule=concentration_of_hco3_rule,
doc="constraint concentration of hco3",
)
def concentration_of_nh3_rule(self):
return (
self.pK_a_IN
== pyo.log10(self.conc_mol_nh4 / (pyo.units.kmole / pyo.units.m**3))
- pyo.log10(self.conc_mol_nh3 / (pyo.units.kmole / pyo.units.m**3))
+ self.pH
)
self.concentration_of_nh3 = pyo.Constraint(
rule=concentration_of_nh3_rule,
doc="constraint concentration of nh3",
)
# TO DO: use correct conversion number
def concentration_of_co2_rule(self):
return (
self.conc_mol_co2
== self.conc_mass_comp_ref["S_IC"] / mw_c - self.conc_mol_hco3
)
self.concentration_of_co2 = pyo.Constraint(
rule=concentration_of_co2_rule,
doc="constraint concentration of co2",
)
def concentration_of_nh4_rule(self):
return (
self.conc_mol_nh4
== self.conc_mass_comp_ref["S_IN"] / mw_n - self.conc_mol_nh3
)
self.concentration_of_nh4 = pyo.Constraint(
rule=concentration_of_nh4_rule,
doc="constraint concentration of nh4",
)
def S_H_rule(self):
return (
self.state_ref.cations
+ self.conc_mol_nh4
+ self.S_H
- self.conc_mol_hco3
- self.conc_mass_ac / (64 * pyo.units.kg / pyo.units.kmol)
- self.conc_mass_pro / (112 * pyo.units.kg / pyo.units.kmol)
- self.conc_mass_bu / (160 * pyo.units.kg / pyo.units.kmol)
- self.conc_mass_va / (208 * pyo.units.kg / pyo.units.kmol)
- 10 ** (self.pH - self.pKW) * (pyo.units.kmole / pyo.units.m**3)
- self.state_ref.anions
== 0
)
self.S_H_cons = pyo.Constraint(
rule=S_H_rule,
doc="constraint concentration of H",
)
def rule_I_IN_lim(self):
return 1 / (
1 + self.params.K_S_IN / (self.conc_mass_comp_ref["S_IN"] / mw_n)
)
self.I_IN_lim = pyo.Expression(
rule=rule_I_IN_lim,
doc="Inhibition function related to secondary substrate; inhibit uptake when inorganic nitrogen S_IN~ 0",
)
def rule_I_h2_fa(self):
return 1 / (1 + self.conc_mass_comp_ref["S_h2"] / self.params.K_I_h2_fa)
self.I_h2_fa = pyo.Expression(
rule=rule_I_h2_fa,
doc="hydrogen inhibition attributed to long chain fatty acids",
)
def rule_I_h2_c4(self):
return 1 / (1 + self.conc_mass_comp_ref["S_h2"] / self.params.K_I_h2_c4)
self.I_h2_c4 = pyo.Expression(
rule=rule_I_h2_c4,
doc="hydrogen inhibition attributed to valerate and butyrate uptake",
)
def rule_I_h2_pro(self):
return 1 / (1 + self.conc_mass_comp_ref["S_h2"] / self.params.K_I_h2_pro)
self.I_h2_pro = pyo.Expression(
rule=rule_I_h2_pro,
doc="hydrogen inhibition attributed to propionate uptake",
)
def rule_I_nh3(self):
return 1 / (1 + self.conc_mol_nh3 / self.params.K_I_nh3)
self.I_nh3 = pyo.Expression(
rule=rule_I_nh3, doc="ammonia inibition attributed to acetate uptake"
)
def rule_I_pH_aa(self):
return (
-3
* (
smooth_max(0, self.params.pH_UL_aa - self.pH, eps=1e-8)
/ (self.params.pH_UL_aa - self.params.pH_LL_aa)
)
** 2
)
self.I_pH_aa = pyo.Expression(
rule=rule_I_pH_aa,
doc="pH inhibition of amino-acid-utilizing microorganisms",
)
def rule_I_pH_ac(self):
return (
-3
* (
smooth_max(0, self.params.pH_UL_ac - self.pH, eps=1e-8)
/ (self.params.pH_UL_ac - self.params.pH_LL_ac)
)
** 2
)
self.I_pH_ac = pyo.Expression(
rule=rule_I_pH_ac, doc="pH inhibition of acetate-utilizing microorganisms"
)
def rule_I_pH_h2(self):
return (
-3
* (
smooth_max(0, self.params.pH_UL_h2 - self.pH, eps=1e-8)
/ (self.params.pH_UL_h2 - self.params.pH_LL_h2)
)
** 2
)
self.I_pH_h2 = pyo.Expression(
rule=rule_I_pH_h2, doc="pH inhibition of hydrogen-utilizing microorganisms"
)
def rule_I(self, r):
if r == "R5" or r == "R6":
return self.I[r] == pyo.exp(self.I_pH_aa) * self.I_IN_lim
elif r == "R7":
return self.I[r] == pyo.exp(self.I_pH_aa) * self.I_IN_lim * self.I_h2_fa
elif r == "R8" or r == "R9":
return self.I[r] == pyo.exp(self.I_pH_aa) * self.I_IN_lim * self.I_h2_c4
elif r == "R10":
return (
self.I[r] == pyo.exp(self.I_pH_aa) * self.I_IN_lim * self.I_h2_pro
)
elif r == "R11":
return self.I[r] == pyo.exp(self.I_pH_ac) * self.I_IN_lim * self.I_nh3
elif r == "R12":
return self.I[r] == pyo.exp(self.I_pH_h2) * self.I_IN_lim
else:
return self.I[r] == 1.0
self.I_fun = pyo.Constraint(
self.params.rate_reaction_idx,
rule=rule_I,
doc="Process inhibition functions",
)
try:
def rate_expression_rule(b, r):
if r == "R1":
# R1: Disintegration
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dis * b.conc_mass_comp_ref["X_c"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R2":
# R2: Hydrolysis of carbohydrates
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_hyd_ch * b.conc_mass_comp_ref["X_ch"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R3":
# R3: Hydrolysis of proteins
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_hyd_pr * b.conc_mass_comp_ref["X_pr"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R4":
# R4: Hydrolysis of lipids
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_hyd_li * b.conc_mass_comp_ref["X_li"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R5":
# R5: Uptake of sugars
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_su
* b.conc_mass_comp_ref["S_su"]
/ (b.params.K_S_su + b.conc_mass_comp_ref["S_su"])
* b.conc_mass_comp_ref["X_su"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R6":
# R6: Uptake of amino acids
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_aa
* b.conc_mass_comp_ref["S_aa"]
/ (b.params.K_S_aa + b.conc_mass_comp_ref["S_aa"])
* b.conc_mass_comp_ref["X_aa"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R7":
# R7: Uptake of long chain fatty acids (LCFAs)
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_fa
* b.conc_mass_comp_ref["S_fa"]
/ (b.params.K_S_fa + b.conc_mass_comp_ref["S_fa"])
* b.conc_mass_comp_ref["X_fa"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R8":
# R8: Uptake of valerate
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_c4
* b.conc_mass_comp_ref["S_va"]
/ (b.params.K_S_c4 + b.conc_mass_comp_ref["S_va"])
* b.conc_mass_comp_ref["X_c4"]
* (
b.conc_mass_comp_ref["S_va"]
/ (
b.conc_mass_comp_ref["S_va"]
+ b.conc_mass_comp_ref["S_bu"]
+ 1e-10 * pyo.units.kg / pyo.units.m**3
)
)
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R9":
# R9: Uptake of butyrate
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_c4
* b.conc_mass_comp_ref["S_bu"]
/ (b.params.K_S_c4 + b.conc_mass_comp_ref["S_bu"])
* b.conc_mass_comp_ref["X_c4"]
* (
b.conc_mass_comp_ref["S_bu"]
/ (
b.conc_mass_comp_ref["S_va"]
+ b.conc_mass_comp_ref["S_bu"]
+ 1e-10 * pyo.units.kg / pyo.units.m**3
)
)
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R10":
# R10: Uptake of propionate
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_pro
* b.conc_mass_comp_ref["S_pro"]
/ (b.params.K_S_pro + b.conc_mass_comp_ref["S_pro"])
* b.conc_mass_comp_ref["X_pro"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R11":
# R11: Uptake of acetate
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_ac
* b.conc_mass_comp_ref["S_ac"]
/ (b.params.K_S_ac + b.conc_mass_comp_ref["S_ac"])
* b.conc_mass_comp_ref["X_ac"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R12":
# R12: Uptake of hydrogen
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_m_h2
* b.conc_mass_comp_ref["S_h2"]
/ (b.params.K_S_h2 + b.conc_mass_comp_ref["S_h2"])
* b.conc_mass_comp_ref["X_h2"]
* b.I[r],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R13":
# R13: Decay of X_su
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_su * b.conc_mass_comp_ref["X_su"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R14":
# R14: Decay of X_aa
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_aa * b.conc_mass_comp_ref["X_aa"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R15":
# R15: Decay of X_fa
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_fa * b.conc_mass_comp_ref["X_fa"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R16":
# R16: Decay of X_c4
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_c4 * b.conc_mass_comp_ref["X_c4"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R17":
# R17: Decay of X_pro
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_pro * b.conc_mass_comp_ref["X_pro"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R18":
# R18: Decay of X_ac
return b.reaction_rate[r] == pyo.units.convert(
b.params.k_dec_X_ac * b.conc_mass_comp_ref["X_ac"],
to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s,
)
elif r == "R19":
# R19: Decay of X_h2
return b.reaction_rate[r] == (
pyo.units.convert(b.params.k_dec_X_h2, to_units=1 / pyo.units.s)
* b.conc_mass_comp_ref["X_h2"]
)
else:
raise BurntToast()
self.rate_expression = pyo.Constraint(
self.params.rate_reaction_idx,
rule=rate_expression_rule,
doc="ADM1 rate expressions",
)
except AttributeError:
# If constraint fails, clean up so that DAE can try again later
self.del_component(self.reaction_rate)
self.del_component(self.rate_expression)
raise
for i, c in self.rates.items():
iscale.set_scaling_factor(self.reaction_rate[i], 1 / c)
iscale.set_scaling_factor(self.I, 1e1)
iscale.set_scaling_factor(self.conc_mass_va, 1e2)
iscale.set_scaling_factor(self.conc_mass_bu, 1e2)
iscale.set_scaling_factor(self.conc_mass_pro, 1e2)
iscale.set_scaling_factor(self.conc_mass_ac, 1e1)
iscale.set_scaling_factor(self.conc_mol_hco3, 1e1)
iscale.set_scaling_factor(self.conc_mol_nh3, 1e3)
iscale.set_scaling_factor(self.conc_mol_co2, 1e3)
iscale.set_scaling_factor(self.conc_mol_nh4, 1e1)
iscale.set_scaling_factor(self.S_H, 1e5)
iscale.set_scaling_factor(self.pKW, 1e0)
iscale.set_scaling_factor(self.pK_a_co2, 1e0)
iscale.set_scaling_factor(self.pK_a_IN, 1e0)
iscale.set_scaling_factor(self.pH, 1e0)
[docs] def get_reaction_rate_basis(self):
return MaterialFlowBasis.mass
def calculate_scaling_factors(self):
super().calculate_scaling_factors()
for i, c in self.rate_expression.items():
iscale.constraint_scaling_transform(
c,
iscale.get_scaling_factor(
self.reaction_rate[i], default=1, warning=True
),
overwrite=True,
)