From 624c3f67adc20ba23f7f69153a7fd96463da6170 Mon Sep 17 00:00:00 2001 From: Yalin Date: Wed, 4 Oct 2023 23:47:31 -0400 Subject: [PATCH 01/16] take interface test out, will use bsm2 when it's done --- tests/test_exposan.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_exposan.py b/tests/test_exposan.py index 50062fcc..b716f8d5 100644 --- a/tests/test_exposan.py +++ b/tests/test_exposan.py @@ -37,15 +37,15 @@ def test_exposan(): bsm1_sys.simulate(t_span=(0,10), method='BDF') print(get_SRT(bsm1_sys, biomass_IDs=biomass_IDs)) # to test the `get_SRT` function + #!!! Will use bsm2 to test the junction models + # from exposan.interface import create_system as create_inter_system + # sys_inter = create_inter_system() + # sys_inter.simulate(method='BDF', t_span=(0, 3)) # the default 'RK45' method can't solve it + from exposan.cas import create_system as create_cas_system cas_sys = create_cas_system() cas_sys.simulate() - from exposan.interface import create_system as create_inter_system - sys_inter = create_inter_system() - #!!! Temporarily disable this test while trying to fix the issue - # sys_inter.simulate(method='BDF', t_span=(0, 3)) # the default 'RK45' method can't solve it - ##### Systems with costs/impacts ##### from qsdsan.utils import clear_lca_registries From 9263bc8d34fd15dc3b75ea35617789ae01a71464 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 5 Oct 2023 08:58:50 -0700 Subject: [PATCH 02/16] modified ADM1 complete component set --- qsdsan/_component.py | 23 ++--- qsdsan/processes/__init__.py | 3 + qsdsan/processes/_madm1.py | 161 +++++++++++++++++++++++++++++++++++ qsdsan/utils/cod.py | 10 ++- 4 files changed, 183 insertions(+), 14 deletions(-) create mode 100644 qsdsan/processes/_madm1.py diff --git a/qsdsan/_component.py b/qsdsan/_component.py index 1a2d60a9..98f43b0a 100644 --- a/qsdsan/_component.py +++ b/qsdsan/_component.py @@ -233,6 +233,10 @@ def __new__(cls, ID, search_ID=None, formula=None, phase=None, measured_as=None, self._chem_MW = molecular_weight(self.atoms) if phase: lock_phase(self, phase) + self._particle_size = particle_size + self._degradability = degradability + self._organic = organic + self.description = description self._measured_as = measured_as self.i_mass = i_mass self.i_C = i_C @@ -245,10 +249,7 @@ def __new__(cls, ID, search_ID=None, formula=None, phase=None, measured_as=None, self.f_BOD5_COD = f_BOD5_COD self.f_uBOD_COD = f_uBOD_COD self.f_Vmass_Totmass = f_Vmass_Totmass - self._particle_size = particle_size - self._degradability = degradability - self._organic = organic - self.description = description + if not self.MW and not self.formula: self.MW = 1. self.i_COD = i_COD self.i_NOD = i_NOD @@ -340,7 +341,8 @@ def i_mass(self, i): chem_charge = charge_from_formula(self.formula) Cr2O7 = - cod_test_stoichiometry(self.atoms, chem_charge)['Cr2O7-2'] cod = Cr2O7 * 1.5 * molecular_weight({'O':2}) - i = self.chem_MW/cod + try: i = self.chem_MW/cod + except: breakpoint() elif self.measured_as: raise AttributeError(f'Must specify i_mass for component {self.ID} ' f'measured as {self.measured_as}.') @@ -532,7 +534,7 @@ def i_COD(self): def i_COD(self, i): if i is not None: self._i_COD = check_return_property('i_COD', i) else: - if self.organic or self.formula in ('H2', 'O2', 'N2', 'NO2-', 'NO3-'): + if self.organic or self.formula in ('H2', 'O2', 'N2', 'NO2-', 'NO3-', 'H2S', 'S'): if self.measured_as == 'COD': self._i_COD = 1. elif not self.atoms: raise AttributeError(f"Must specify `i_COD` for organic component {self.ID}, " @@ -792,6 +794,10 @@ def from_chemical(cls, ID, chemical=None, formula=None, phase=None, measured_as= new.Tm, new.Tb, new.eos, new.phase_ref, new.S0) TDependentProperty.RAISE_PROPERTY_CALCULATION_ERROR = True + new.description = description + new._particle_size = particle_size + new._degradability = degradability + new._organic = organic new._measured_as = measured_as new.i_mass = i_mass new.i_C = i_C @@ -804,10 +810,7 @@ def from_chemical(cls, ID, chemical=None, formula=None, phase=None, measured_as= new.f_BOD5_COD = f_BOD5_COD new.f_uBOD_COD = f_uBOD_COD new.f_Vmass_Totmass = f_Vmass_Totmass - new.description = description - new._particle_size = particle_size - new._degradability = degradability - new._organic = organic + new.i_COD = i_COD new.i_NOD = i_NOD return new \ No newline at end of file diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index a8f62825..9e0c5e9d 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -17,6 +17,7 @@ from ._asm1 import * from ._asm2d import * from ._adm1 import * +from ._madm1 import * from ._decay import * from ._kinetic_reaction import * from ._pm2 import * @@ -26,6 +27,7 @@ _asm1, _asm2d, _adm1, + _madm1, _decay, _kinetic_reaction, ) @@ -35,6 +37,7 @@ *_asm1.__all__, *_asm2d.__all__, *_adm1.__all__, + *_madm1.__all__, *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py new file mode 100644 index 00000000..6a819f4f --- /dev/null +++ b/qsdsan/processes/_madm1.py @@ -0,0 +1,161 @@ +# -*- coding: utf-8 -*- +''' +QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems + +This module is developed by: + Joy Zhang + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt +for license details. +''' + +# from thermosteam.utils import chemicals_user +from chemicals.elements import molecular_weight as get_mw +from qsdsan import Component, Components, Process, Processes, CompiledProcesses +import numpy as np, qsdsan.processes as pc, qsdsan as qs +# from qsdsan.utils import ospath, data_path +# from scipy.optimize import brenth +# from warnings import warn + +# from qsdsan.processes import create_adm1_cmps + + +__all__ = ('create_madm1_cmps', + ) + +#%% +# C_mw = get_mw({'C':1}) +# N_mw = get_mw({'N':1}) +# P_mw = get_mw({'P':1}) +# S_mw = get_mw({'S':1}) +Fe_mw = get_mw({'Fe':1}) +O_mw = get_mw({'O':1}) + +def create_madm1_cmps(set_thermo=True): + + # Components from the original ADM1 + # ********************************* + _cmps = pc.create_adm1_cmps(False) + S_aa = _cmps.S_aa + X_pr = _cmps.X_pr + S_aa.i_C = X_pr.i_C = 0.36890 + S_aa.i_N = X_pr.i_N = 0.11065 + S_aa.i_mass = X_pr.i_mass = 1/1.35566 + + S_fa = _cmps.S_fa + S_fa.formula = 'C25H52O3' + S_bu = _cmps.S_bu + S_bu.formula = 'C4H8O2' + S_pro = _cmps.S_pro + S_pro.formula = 'C3H6O2' + S_ac = _cmps.S_ac + S_ac.formula = 'C2H4O2' + + S_I = _cmps.S_I + X_I = _cmps.X_I + S_I.i_C = X_I.i_C = 0.36178 + S_I.i_N = X_I.i_N = 0.06003 + S_I.i_P = X_I.i_P = 0.00649 + S_I.i_mass = X_I.i_mass = 1/1.54100 + + X_ch = _cmps.X_ch + X_ch.formula = 'C24H48O24' + # _cmps.X_li.formula = 'C64H119O7.5P' + X_li = X_pr.copy('X_li') + X_li.i_C = 0.26311 + X_li.i_N = 0. + X_li.i_P = 0.01067 + X_li.i_mass = 1/2.81254 + + adm1_biomass = (_cmps.X_su, _cmps.X_aa, _cmps.X_fa, _cmps.X_c4, _cmps.X_pro, _cmps.X_ac, _cmps.X_h2) + for bio in adm1_biomass: + # bio.formula = 'C5H7O2NP0.113' + bio.i_C = 0.36612 + bio.i_N = 0.08615 + bio.i_P = 0.02154 + bio.i_mass = 1/1.39300 + + # P related components from ASM2d + # ******************************* + asm_cmps = pc.create_asm2d_cmps(False) + X_PHA = asm_cmps.X_PHA + X_PHA.formula = '(C2H4O)n' + # X_PHA.i_C = 0.3 + # X_PHA.i_mass = 0.55 + + X_PAO = _cmps.X_su.copy('X_PAO') + X_PAO.description = 'Phosphorus-accumulating organism biomass' + + # Additional components for P, S, Fe extensions + # ********************************************* + S_IP = asm_cmps.S_PO4.copy('S_IP') + S_K = Component.from_chemical('S_K', chemical='K+', + description='Potassium ion', + measured_as='K', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_Mg = Component.from_chemical('S_Mg', chemical='Mg2+', + description='Magnesium ion', + measured_as='Mg', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_SO4 = Component.from_chemical('S_SO4', chemical='SO4-2', + description='Sulfate', + measured_as='S', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_IS = Component.from_chemical('S_IS', chemical='H2S', + description='Hydrogen sulfide', + measured_as='COD', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + + X_hSRB = X_PAO.copy('X_hSRB') + X_hSRB.description = 'sulfate-reducing biomass, utilizing H2' + X_aSRB = X_PAO.copy('X_aSRB') + X_aSRB.description = 'sulfate-reducing biomass, utilizing acetate' + X_pSRB = X_PAO.copy('X_pSRB') + X_pSRB.description = 'sulfate-reducing biomass, utilizing propionate' + X_c4SRB = X_PAO.copy('X_c4SRB') + X_c4SRB.description = 'sulfate-reducing biomass, utilizing butyrate and valerate' + + S_S0 = Component.from_chemical('S_S0', chemical='S', + description='Sulfur', + measured_as='COD', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_Fe3 = Component.from_chemical('S_Fe3', chemical='Fe3+', + description='Iron (III)', + measured_as='Fe', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_Fe2 = Component.from_chemical('S_Fe2', chemical='Fe2+', + description='Iron (II)', + measured_as='Fe', + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_Fe2.i_COD = 0.5*O_mw/Fe_mw + S_Fe2.measured_as = 'COD' + + cmps_madm1 = Components([_cmps.S_su, S_aa, S_fa, _cmps.S_va, S_bu, + S_pro, S_ac, _cmps.S_h2, _cmps.S_ch4, + _cmps.S_IC, _cmps.S_IN, S_IP, S_I, + X_ch, X_pr, X_li, *adm1_biomass, X_I, + X_PHA, asm_cmps.X_PP, X_PAO, S_K, S_Mg, + S_SO4, S_IS, X_hSRB, X_aSRB, X_pSRB, X_c4SRB, + S_S0, S_Fe3, S_Fe2, + _cmps.H2O]) + cmps_madm1.default_compile() + + if set_thermo: qs.set_thermo(cmps_madm1) + return cmps_madm1 + +#%% diff --git a/qsdsan/utils/cod.py b/qsdsan/utils/cod.py index 7d4b260b..965bbe99 100644 --- a/qsdsan/utils/cod.py +++ b/qsdsan/utils/cod.py @@ -107,10 +107,12 @@ def cod_test_stoichiometry(atoms, charge=0, MW=None, missing_handling='elemental nC, nH, nO, nN, nS, nP = get_CHONSP(atoms) ne = - charge - if nC <= 0 or nH <= 0: - if not (len(atoms) == 1 and nH == 2): # H2 - return {'Cr2O7-2': 0.} - + # if nC <= 0 or nH <= 0: + # if not (len(atoms) == 1 and nH == 2): # H2 + # return {'Cr2O7-2': 0.} + if nC + nH + nS + nP <= 0: + return {'Cr2O7-2': 0.} + nCO2 = nC nNH4 = nN nSO4 = nS From 7c45160d1f780bdc087ce6cf41573122ac115fa4 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 5 Oct 2023 14:51:39 -0700 Subject: [PATCH 03/16] modified ADM1 stoichiometry --- qsdsan/data/process_data/_madm1.tsv | 37 +++++++++ qsdsan/processes/_madm1.py | 122 ++++++++++++++++++++++++++-- 2 files changed, 152 insertions(+), 7 deletions(-) create mode 100644 qsdsan/data/process_data/_madm1.tsv diff --git a/qsdsan/data/process_data/_madm1.tsv b/qsdsan/data/process_data/_madm1.tsv new file mode 100644 index 00000000..2e124782 --- /dev/null +++ b/qsdsan/data/process_data/_madm1.tsv @@ -0,0 +1,37 @@ + S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg S_SO4 S_IS X_hSRB X_aSRB X_pSRB X_c4SRB S_S0 S_Fe3 S_Fe2 +hydrolysis_carbs 1 ? ? ? -1 +hydrolysis_proteins 1 ? ? ? -1 +hydrolysis_lipids 1-f_fa_li f_fa_li ? ? ? -1 +uptake_sugars -1 (1-Y_su)*f_bu_su (1-Y_su)*f_pro_su (1-Y_su)*f_ac_su (1-Y_su)*f_h2_su ? ? ? Y_su +uptake_amino_acids -1 (1-Y_aa)*f_va_aa (1-Y_aa)*f_bu_aa (1-Y_aa)*f_pro_aa (1-Y_aa)*f_ac_aa (1-Y_aa)*f_h2_aa ? ? ? Y_aa +uptake_LCFA -1 (1-Y_fa)*f_ac_fa (1-Y_fa)*f_h2_fa ? ? ? Y_fa +uptake_valerate -1 (1-Y_c4)*f_pro_va (1-Y_c4)*f_ac_va (1-Y_c4)*f_h2_va ? ? ? Y_c4 +uptake_butyrate -1 (1-Y_c4)*f_ac_bu (1-Y_c4)*f_h2_bu ? ? ? Y_c4 +uptake_propionate -1 (1-Y_pro)*f_ac_pro (1-Y_pro)*f_h2_pro ? ? ? Y_pro +uptake_acetate -1 1-Y_ac ? ? ? Y_ac +uptake_h2 -1 1-Y_h2 ? ? ? Y_h2 +decay_Xsu ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xaa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xfa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xc4 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xpro ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xac ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xh2 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +storage_PHA_valerate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_butyrate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_propionate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_actate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +lysis_XPAO ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +lysis_XPP ? ? ? -1 K_XPP Mg_XPP +lysis_XPHA f_va_pha f_bu_pha f_pro_pha f_ac_pha ? ? ? -1 +growth_SRB_h2 -1 ? ? ? -1*(1-Y_hSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_hSRB Y_hSRB +decay_XhSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_acetate -1 ? ? ? -1*(1-Y_aSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_aSRB Y_aSRB +decay_XaSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_propionate -1 (1-Y_pSRB)*f_ac_pro ? ? ? -1*(1-Y_pSRB)*f_is_pro*i_mass_IS*(MW_S0/MW_IS) (1-Y_pSRB)*f_is_pro Y_pSRB +decay_XpSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_butyrate -1 (1-Y_c4SRB)*f_ac_bu ? ? ? -1*(1-Y_c4SRB)*f_is_bu*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_bu Y_c4SRB +growth_SRB_valerate -1 (1-Y_c4SRB)*f_pro_va (1-Y_c4SRB)*f_ac_va ? ? ? -1*(1-Y_c4SRB)*f_is_va*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_va Y_c4SRB +decay_Xc4SRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +Fe3_reduction_h2 -1 -1*i_mass_Fe2 1 +Fe3_reduction_IS -1 (MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1) 1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 6a819f4f..0f991cc9 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -10,21 +10,22 @@ for license details. ''' -# from thermosteam.utils import chemicals_user +from thermosteam.utils import chemicals_user +from thermosteam import settings from chemicals.elements import molecular_weight as get_mw from qsdsan import Component, Components, Process, Processes, CompiledProcesses import numpy as np, qsdsan.processes as pc, qsdsan as qs -# from qsdsan.utils import ospath, data_path +from qsdsan.utils import ospath, data_path # from scipy.optimize import brenth # from warnings import warn -# from qsdsan.processes import create_adm1_cmps +__all__ = ('create_madm1_cmps', 'ModifiedADM1') -__all__ = ('create_madm1_cmps', - ) +_path = ospath.join(data_path, 'process_data/_madm1.tsv') -#%% + +#%% components # C_mw = get_mw({'C':1}) # N_mw = get_mw({'N':1}) # P_mw = get_mw({'P':1}) @@ -41,6 +42,7 @@ def create_madm1_cmps(set_thermo=True): X_pr = _cmps.X_pr S_aa.i_C = X_pr.i_C = 0.36890 S_aa.i_N = X_pr.i_N = 0.11065 + S_aa.i_P = X_pr.i_P = 0. S_aa.i_mass = X_pr.i_mass = 1/1.35566 S_fa = _cmps.S_fa @@ -158,4 +160,110 @@ def create_madm1_cmps(set_thermo=True): if set_thermo: qs.set_thermo(cmps_madm1) return cmps_madm1 -#%% +#%% rate functions + + +#%% modified ADM1 class +_load_components = settings.get_default_chemicals + +@chemicals_user +class ModifiedADM1(CompiledProcesses): + + _cmp_dependent_stoichio = ('K_XPP', 'Mg_XPP', + 'MW_S0', 'MW_IS', + 'i_mass_S0', 'i_mass_IS', 'i_mass_Fe2') + _stoichio_params = (*pc.ADM1._stoichio_params[5:], + 'f_ch_xb', 'f_pr_xb', 'f_li_xb', 'f_xI_xb', 'f_sI_xb', + 'f_va_pha', 'f_bu_pha', 'f_pro_pha', 'f_ac_pha', + 'f_is_pro', 'f_is_bu', 'f_is_va', + 'Y_PO4', 'Y_hSRB', 'Y_aSRB', 'Y_pSRB', 'Y_c4SRB', + *_cmp_dependent_stoichio + ) + # _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'pH_ULs', 'pH_LLs', + # 'KS_IN', 'KI_nh3', 'KIs_h2', + # 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', + # 'T_base', 'components', 'root') + # _acid_base_pairs = (('H+', 'OH-'), ('NH4+', 'NH3'), ('CO2', 'HCO3-'), + # ('HAc', 'Ac-'), ('HPr', 'Pr-'), + # ('HBu', 'Bu-'), ('HVa', 'Va-')) + _biogas_IDs = (*pc.ADM1._biogas_IDs, 'S_IS') + _biomass_IDs = (*pc.ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') + + def __new__(cls, components=None, path=None, + f_ch_xb=0.275, f_pr_xb=0.275, f_li_xb=0.35, f_xI_xb=0.1, + f_fa_li=0.95, f_bu_su=0.1328, f_pro_su=0.2691, f_ac_su=0.4076, + f_va_aa=0.23, f_bu_aa=0.26, f_pro_aa=0.05, f_ac_aa=0.4, + f_ac_fa=0.7, f_pro_va=0.54, f_ac_va=0.31, f_ac_bu=0.8, f_ac_pro=0.57, + Y_su=0.1, Y_aa=0.08, Y_fa=0.06, Y_c4=0.06, Y_pro=0.04, Y_ac=0.05, Y_h2=0.06, + f_va_pha=0.1, f_bu_pha=0.1, f_pro_pha=0.4, + Y_PO4=0.4, Y_hSRB=0.05, Y_aSRB=0.05, Y_pSRB=0.04, Y_c4SRB=0.06, + + q_dis=0.5, q_ch_hyd=10, q_pr_hyd=10, q_li_hyd=10, + k_su=30, k_aa=50, k_fa=6, k_c4=20, k_pro=13, k_ac=8, k_h2=35, + K_su=0.5, K_aa=0.3, K_fa=0.4, K_c4=0.2, K_pro=0.1, K_ac=0.15, K_h2=7e-6, + b_su=0.02, b_aa=0.02, b_fa=0.02, b_c4=0.02, b_pro=0.02, b_ac=0.02, b_h2=0.02, + KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, + pH_limits_aa=(4,5.5), pH_limits_ac=(6,7), pH_limits_h2=(5,6), + T_base=298.15, pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], + Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0], + kLa=200, K_H_base=[7.8e-4, 1.4e-3, 3.5e-2], + K_H_dH=[-4180, -14240, -19410], + **kwargs): + + cmps = _load_components(components) + + if not path: path = _path + self = Processes.load_from_file(path, + components=cmps, + conserved_for=('C', 'N', 'P'), + # parameters=_stoichio_params, + parameters=cls._stoichio_params, + compile=False) + + gas_transfer = [] + # for i in _biogas_IDs: + for i in cls._biogas_IDs: + new_p = Process('%s_transfer' % i.lstrip('S_'), + reaction={i:-1}, + ref_component=i, + conserved_for=(), + parameters=()) + gas_transfer.append(new_p) + self.extend(gas_transfer) + self.compile(to_class=cls) + + stoichio_vals = (f_fa_li, f_bu_su, f_pro_su, f_ac_su, 1-f_bu_su-f_pro_su-f_ac_su, + f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa, 1-f_va_aa-f_bu_aa-f_pro_aa-f_ac_aa, + f_ac_fa, 1-f_ac_fa, f_pro_va, f_ac_va, 1-f_pro_va-f_ac_va, + f_ac_bu, 1-f_ac_bu, f_ac_pro, 1-f_ac_pro, + Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2, + f_ch_xb, f_pr_xb, f_li_xb, f_xI_xb, round(1.0-f_ch_xb-f_pr_xb-f_li_xb-f_xI_xb, 4), + f_va_pha, f_bu_pha, f_pro_pha, 1-f_va_pha-f_bu_pha-f_pro_pha, + 1-f_ac_pro, 1-f_ac_bu, 1-f_pro_va-f_ac_va, + Y_PO4, Y_hSRB, Y_aSRB, Y_pSRB, Y_c4SRB, + cmps.X_PP.i_K, cmps.X_PP.i_Mg, + cmps.S_S0.chem_MW, cmps.S_IS.chem_MW, + cmps.S_S0.i_mass, cmps.S_IS.i_mass, cmps.S_Fe2.i_mass) + # pH_LLs = np.array([pH_limits_aa[0]]*6 + [pH_limits_ac[0], pH_limits_h2[0]]) + # pH_ULs = np.array([pH_limits_aa[1]]*6 + [pH_limits_ac[1], pH_limits_h2[1]]) + # ks = np.array((q_dis, q_ch_hyd, q_pr_hyd, q_li_hyd, + # k_su, k_aa, k_fa, k_c4, k_c4, k_pro, k_ac, k_h2, + # b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2)) + # Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2)) + # KIs_h2 = np.array((KI_h2_fa, KI_h2_c4, KI_h2_c4, KI_h2_pro)) + # K_H_base = np.array(K_H_base) + # K_H_dH = np.array(K_H_dH) + # Ka_base = np.array([10**(-pKa) for pKa in pKa_base]) + # Ka_dH = np.array(Ka_dH) + # root = TempState() + dct = self.__dict__ + dct.update(kwargs) + + # self.set_rate_function(rhos_adm1) + dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) + # self.rate_function._params = dict(zip(cls._kinetic_params, + # [ks, Ks, pH_ULs, pH_LLs, KS_IN*N_mw, + # KI_nh3, KIs_h2, Ka_base, Ka_dH, + # K_H_base, K_H_dH, kLa, + # T_base, self._components, root])) + return self \ No newline at end of file From f72ae9496018bfa7882427c9ea0f77cfec186896 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 5 Oct 2023 18:27:03 -0700 Subject: [PATCH 04/16] modified ADM1 kinetic parameters --- qsdsan/processes/_madm1.py | 98 ++++++++++++++++++++++++-------------- 1 file changed, 61 insertions(+), 37 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 0f991cc9..acf5e913 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -27,9 +27,9 @@ #%% components # C_mw = get_mw({'C':1}) -# N_mw = get_mw({'N':1}) -# P_mw = get_mw({'P':1}) -# S_mw = get_mw({'S':1}) +N_mw = get_mw({'N':1}) +P_mw = get_mw({'P':1}) +S_mw = get_mw({'S':1}) Fe_mw = get_mw({'Fe':1}) O_mw = get_mw({'O':1}) @@ -162,6 +162,8 @@ def create_madm1_cmps(set_thermo=True): #%% rate functions +def rhos_madm1(): + pass #%% modified ADM1 class _load_components = settings.get_default_chemicals @@ -179,15 +181,17 @@ class ModifiedADM1(CompiledProcesses): 'Y_PO4', 'Y_hSRB', 'Y_aSRB', 'Y_pSRB', 'Y_c4SRB', *_cmp_dependent_stoichio ) - # _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'pH_ULs', 'pH_LLs', - # 'KS_IN', 'KI_nh3', 'KIs_h2', - # 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', - # 'T_base', 'components', 'root') - # _acid_base_pairs = (('H+', 'OH-'), ('NH4+', 'NH3'), ('CO2', 'HCO3-'), - # ('HAc', 'Ac-'), ('HPr', 'Pr-'), - # ('HBu', 'Bu-'), ('HVa', 'Va-')) + _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'K_PP', 'K_so4', + 'pH_limits', 'KS_IN', 'KI_nh3', 'KIs_h2', + 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', + 'T_base', 'components', 'root', + ) + _acid_base_pairs = pc.ADM1._acid_base_pairs _biogas_IDs = (*pc.ADM1._biogas_IDs, 'S_IS') _biomass_IDs = (*pc.ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') + _T_base = 298.15 + _K_H_base = [7.8e-4, 1.4e-3, 3.5e-2, 0.105] # biogas species Henry's Law constant [M/bar] + _K_H_dH = [-4180, -14240, -19410, -19180] # Heat of reaction of liquid-gas transfer of biogas species [J/mol] def __new__(cls, components=None, path=None, f_ch_xb=0.275, f_pr_xb=0.275, f_li_xb=0.35, f_xI_xb=0.1, @@ -196,18 +200,25 @@ def __new__(cls, components=None, path=None, f_ac_fa=0.7, f_pro_va=0.54, f_ac_va=0.31, f_ac_bu=0.8, f_ac_pro=0.57, Y_su=0.1, Y_aa=0.08, Y_fa=0.06, Y_c4=0.06, Y_pro=0.04, Y_ac=0.05, Y_h2=0.06, f_va_pha=0.1, f_bu_pha=0.1, f_pro_pha=0.4, - Y_PO4=0.4, Y_hSRB=0.05, Y_aSRB=0.05, Y_pSRB=0.04, Y_c4SRB=0.06, - - q_dis=0.5, q_ch_hyd=10, q_pr_hyd=10, q_li_hyd=10, + Y_PO4=0.4, Y_hSRB=0.05, Y_aSRB=0.05, Y_pSRB=0.04, Y_c4SRB=0.06, + q_ch_hyd=10, q_pr_hyd=10, q_li_hyd=10, k_su=30, k_aa=50, k_fa=6, k_c4=20, k_pro=13, k_ac=8, k_h2=35, K_su=0.5, K_aa=0.3, K_fa=0.4, K_c4=0.2, K_pro=0.1, K_ac=0.15, K_h2=7e-6, b_su=0.02, b_aa=0.02, b_fa=0.02, b_c4=0.02, b_pro=0.02, b_ac=0.02, b_h2=0.02, - KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, + q_pha=3.0, b_pao=0.2, b_pp=0.2, b_pha=0.2, K_A=4e-3, K_PP=0.01, + k_hSRB=41.125, k_aSRB=10., k_pSRB=16.25, k_c4SRB=23, + b_hSRB=0.02, b_aSRB=0.02, b_pSRB=0.02, b_c4SRB=0.02, + K_hSRB=5.96e-6, K_aSRB=0.176, K_pSRB=0.088, K_c4SRB=0.1739, + K_so4_hSRB=1.04e-4*S_mw, K_so4_aSRB=2e-4*S_mw, K_so4_pSRB=2e-4*S_mw, K_so4_c4SRB=2e-4*S_mw, + k_Fe3t2=1e9/Fe_mw, + KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, KS_IP=2e-5, + KI_h2s_c4=0.481, KI_h2s_pro=0.481, KI_h2s_ac=0.460, KI_h2s_h2=0.400, + KI_h2s_c4SRB=0.520, KI_h2s_pSRB=0.520, KI_h2s_aSRB=0.499, KI_h2s_hSRB=0.499, pH_limits_aa=(4,5.5), pH_limits_ac=(6,7), pH_limits_h2=(5,6), - T_base=298.15, pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], + pH_limits_aa_SRB=(6,7), pH_limits_ac_SRB=(6,7), pH_limits_h2_SRB=(5,6), + kLa=200, + pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0], - kLa=200, K_H_base=[7.8e-4, 1.4e-3, 3.5e-2], - K_H_dH=[-4180, -14240, -19410], **kwargs): cmps = _load_components(components) @@ -216,12 +227,10 @@ def __new__(cls, components=None, path=None, self = Processes.load_from_file(path, components=cmps, conserved_for=('C', 'N', 'P'), - # parameters=_stoichio_params, parameters=cls._stoichio_params, compile=False) gas_transfer = [] - # for i in _biogas_IDs: for i in cls._biogas_IDs: new_p = Process('%s_transfer' % i.lstrip('S_'), reaction={i:-1}, @@ -244,26 +253,41 @@ def __new__(cls, components=None, path=None, cmps.X_PP.i_K, cmps.X_PP.i_Mg, cmps.S_S0.chem_MW, cmps.S_IS.chem_MW, cmps.S_S0.i_mass, cmps.S_IS.i_mass, cmps.S_Fe2.i_mass) - # pH_LLs = np.array([pH_limits_aa[0]]*6 + [pH_limits_ac[0], pH_limits_h2[0]]) - # pH_ULs = np.array([pH_limits_aa[1]]*6 + [pH_limits_ac[1], pH_limits_h2[1]]) - # ks = np.array((q_dis, q_ch_hyd, q_pr_hyd, q_li_hyd, - # k_su, k_aa, k_fa, k_c4, k_c4, k_pro, k_ac, k_h2, - # b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2)) - # Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2)) - # KIs_h2 = np.array((KI_h2_fa, KI_h2_c4, KI_h2_c4, KI_h2_pro)) - # K_H_base = np.array(K_H_base) - # K_H_dH = np.array(K_H_dH) - # Ka_base = np.array([10**(-pKa) for pKa in pKa_base]) - # Ka_dH = np.array(Ka_dH) + + pH_limits = np.array([pH_limits_aa, pH_limits_ac, pH_limits_h2, + pH_limits_h2_SRB, pH_limits_ac_SRB, pH_limits_aa_SRB]).T + + ks = np.array((q_ch_hyd, q_pr_hyd, q_li_hyd, + k_su, k_aa, k_fa, k_c4, k_c4, k_pro, k_ac, k_h2, + b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2, # original ADM1 + q_pha, q_pha, q_pha, q_pha, b_pao, b_pp, b_pha, # P extension + k_hSRB, b_hSRB, k_aSRB, b_aSRB, k_pSRB, b_pSRB, k_c4SRB, k_c4SRB, b_c4SRB, # S extension + k_Fe3t2, k_Fe3t2)) # Fe extension + + Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2, # original ADM1 + K_A, # P extension + K_hSRB, K_aSRB, K_pSRB, K_c4SRB)) # S extension + K_so4 = np.array((K_so4_hSRB, K_so4_aSRB, K_so4_pSRB, K_so4_c4SRB)) + + KIs_h2 = np.array((KI_h2_fa, KI_h2_c4, KI_h2_c4, KI_h2_pro)) + KIs_h2s = np.array((KI_h2s_c4, KI_h2s_c4, KI_h2s_pro, KI_h2s_ac, KI_h2s_h2, + KI_h2s_hSRB, KI_h2s_aSRB, KI_h2s_pSRB, KI_h2s_c4SRB, KI_h2s_c4SRB)) + K_H_base = np.array(cls._K_H_base) + K_H_dH = np.array(cls._K_H_dH) + Ka_base = np.array([10**(-pKa) for pKa in pKa_base]) + Ka_dH = np.array(Ka_dH) # root = TempState() dct = self.__dict__ dct.update(kwargs) - - # self.set_rate_function(rhos_adm1) + + self.set_rate_function(rhos_madm1) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) - # self.rate_function._params = dict(zip(cls._kinetic_params, - # [ks, Ks, pH_ULs, pH_LLs, KS_IN*N_mw, - # KI_nh3, KIs_h2, Ka_base, Ka_dH, - # K_H_base, K_H_dH, kLa, - # T_base, self._components, root])) + self.rate_function._params = dict(zip(cls._kinetic_params, + [ks, Ks, K_PP, K_so4, + pH_limits, KS_IN*N_mw, KS_IP*P_mw, + KI_nh3, KIs_h2, KIs_h2s, + Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, + cls.T_base, self._components, + # root, + ])) return self \ No newline at end of file From e2240bff82378365d4f682edaa3f0abb9f4a6e45 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Fri, 6 Oct 2023 15:05:49 -0700 Subject: [PATCH 05/16] add docstring to `ModifiedADM1` --- qsdsan/processes/_madm1.py | 157 +++++++++++++++++++++++++++++++++++-- 1 file changed, 152 insertions(+), 5 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index acf5e913..9d3218d3 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -165,12 +165,161 @@ def create_madm1_cmps(set_thermo=True): def rhos_madm1(): pass + #%% modified ADM1 class _load_components = settings.get_default_chemicals @chemicals_user class ModifiedADM1(CompiledProcesses): - + """ + Modified Anaerobic Digestion Model no.1 [1]_ + + Parameters + ---------- + f_ch_xb : float, optional + Fraction of carbohydrates as biomass decay product. The default is 0.275. + f_pr_xb : flaot, optional + Fraction of proteins as biomass decay product. The default is 0.275. + f_li_xb : float, optional + Fraction of lipids as biomass decay product. The default is 0.35. + f_xI_xb : float, optional + Fraction of inert particulates as biomass decay product. The default is 0.1. + f_va_pha : float, optional + Fraction of valerate as PHA lysis product. The default is 0.1. + f_bu_pha : float, optional + Fraction of butyrate as PHA lysis product. The default is 0.1. + f_pro_pha : float, optional + Fraction of propionate as PHA lysis product. The default is 0.4. + Y_PO4 : float, optional + Poly-phosphorus (PP) required for PHA storage [kg P/kg COD]. The default is 0.4. + Y_hSRB : float, optional + Sulfide-reducing biomass (SRB) yield of hydrogen uptake [kg COD/kg COD]. + The default is 0.05. + Y_aSRB : float, optional + SRB yield of acetate uptake [kg COD/kg COD]. The default is 0.05. + Y_pSRB : float, optional + SRB yield of propionate uptake [kg COD/kg COD]. The default is 0.04. + Y_c4SRB : float, optional + SRB yield of butyrate or valerate uptake [kg COD/kg COD]. + The default is 0.06. + q_pha : float, optional + Maximum specific rate constant for PHA storage by phosphorus-accumulating + organisms (PAOs) [d^(-1)]. The default is 3.0. + b_pao : float, optional + PAO lysis rate constant [d^(-1)]. The default is 0.2. + b_pp : float, optional + PP lysis rate constant [d^(-1)]. The default is 0.2. + b_pha : float, optional + PHA lysis rate constant [d^(-1)]. The default is 0.2. + K_A : float, optional + Substrate half saturation coefficient for PHA storage [kg COD/m3]. + The default is 4e-3. + K_PP : float, optional + PP half saturation coefficient for PHA storage [kg P (X_PP)/kg COD (X_PHA)]. + The default is 0.01. + k_hSRB : float, optional + Maximum specific growth rate constant of hydrogen-uptaking SRB [d^(-1)]. + The default is 41.125. + k_aSRB : float, optional + Maximum specific growth rate constant of acetate-uptaking SRB [d^(-1)]. + The default is 10.. + k_pSRB : float, optional + Maximum specific growth rate constant of propionate-uptaking SRB [d^(-1)]. + The default is 16.25. + k_c4SRB : float, optional + Maximum specific growth rate constant of butyrate- or valerate-uptaking + SRB [d^(-1)]. The default is 23. + b_hSRB : float, optional + Hydrogen-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. + b_aSRB : float, optional + Acetate-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. + b_pSRB : float, optional + Propionate-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. + b_c4SRB : float, optional + Butyrate- or valerate-uptaking SRB decay rate constant [d^(-1)]. + The default is 0.02. + K_hSRB : float, optional + Substrate half saturation coefficient of hydrogen uptake by SRB + [kg COD/m3]. The default is 5.96e-6. + K_aSRB : float, optional + Substrate half saturation coefficient of acetate uptake by SRB + [kg COD/m3]. The default is 0.176. + K_pSRB : float, optional + Substrate half saturation coefficient of propionate uptake by SRB + [kg COD/m3]. The default is 0.088. + K_c4SRB : float, optional + Substrate half saturation coefficient of butyrate or valerate uptake by + SRB [kg COD/m3]. The default is 0.1739. + K_so4_hSRB : float, optional + Sulfate half saturation coefficient of SRB uptaking hydrogen [kg S/m3]. + The default is 3.335e-3. + K_so4_aSRB : float, optional + Sulfate half saturation coefficient of SRB uptaking acetate [kg S/m3]. + The default is 6.413e-3. + K_so4_pSRB : float, optional + Sulfate half saturation coefficient of SRB uptaking propionate [kg S/m3]. + The default is 6.413e-3. + K_so4_c4SRB : float, optional + Sulfate half saturation coefficient of SRB uptaking butyrate or valerate + [kg S/m3]. The default is 6.413e-3. + k_Fe3t2 : float, optional + Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)]. + The default is 1.79e7. + KS_IP : float, optional + Inorganic phosphorus (nutrient) inhibition coefficient for soluble + substrate uptake [M]. The default is 2e-5. + KI_h2s_c4 : float, optional + H2S half inhibition coefficient for butyrate or valerate uptake + [kg COD/m3]. The default is 0.481. + KI_h2s_pro : float, optional + H2S half inhibition coefficient for propionate uptake [kg COD/m3]. + The default is 0.481. + KI_h2s_ac : float, optional + H2S half inhibition coefficient for acetate uptake [kg COD/m3]. + The default is 0.460. + KI_h2s_h2 : float, optional + H2S half inhibition coefficient for hydrogen uptake [kg COD/m3]. + The default is 0.400. + KI_h2s_c4SRB : float, optional + H2S half inhibition coefficient for butyrate or valerate uptake by SRB + [kg COD/m3]. The default is 0.520. + KI_h2s_pSRB : float, optional + H2S half inhibition coefficient for propionate uptake by SRB [kg COD/m3]. + The default is 0.520. + KI_h2s_aSRB : float, optional + H2S half inhibition coefficient for acetate uptake by SRB [kg COD/m3]. + The default is 0.499. + KI_h2s_hSRB : float, optional + H2S half inhibition coefficient for hydrogen uptake by SRB [kg COD/m3]. + The default is 0.499. + pH_limits_aa_SRB : 2-tuple, optional + Lower and upper limits of pH inhibition for acetogenosis by SRB, + unitless. The default is (6,7). + pH_limits_ac_SRB : 2-tuple, optional + Lower and upper limits of pH inhibition for acetate uptake by SRB, + unitless. The default is (6,7). + pH_limits_h2_SRB : 2-tuple, optional + Lower and upper limits of pH inhibition for hydrogen uptake by SRB, + unitless. The default is (5,6). + + Examples + -------- + ... + + References + ---------- + .. [1] Flores-Alsina, X., Solon, K., Kazadi Mbamba, C., Tait, S., + Gernaey, K. V., Jeppsson, U., ; Batstone, D. J. (2016). + Modelling phosphorus (P), sulfur (S) and iron (Fe) interactions + for dynamic simulations of anaerobic digestion processes. + Water Research, 95, 370–382. https://doi.org/10.1016/J.WATRES.2016.03.012 + + See Also + -------- + `qsdsan.processes.ADM1 `_ + + """ + _cmp_dependent_stoichio = ('K_XPP', 'Mg_XPP', 'MW_S0', 'MW_IS', 'i_mass_S0', 'i_mass_IS', 'i_mass_Fe2') @@ -216,10 +365,8 @@ def __new__(cls, components=None, path=None, KI_h2s_c4SRB=0.520, KI_h2s_pSRB=0.520, KI_h2s_aSRB=0.499, KI_h2s_hSRB=0.499, pH_limits_aa=(4,5.5), pH_limits_ac=(6,7), pH_limits_h2=(5,6), pH_limits_aa_SRB=(6,7), pH_limits_ac_SRB=(6,7), pH_limits_h2_SRB=(5,6), - kLa=200, - pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], - Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0], - **kwargs): + kLa=200, pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], + Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0],**kwargs): cmps = _load_components(components) From f226608916c6ab882d3e7e70ba982930689c4c31 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Fri, 6 Oct 2023 17:12:07 -0700 Subject: [PATCH 06/16] mADM1 kinetic rate equations --- qsdsan/processes/_madm1.py | 156 ++++++++++++++++++++++++++++++++++--- 1 file changed, 147 insertions(+), 9 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 9d3218d3..b8f83db6 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -16,6 +16,16 @@ from qsdsan import Component, Components, Process, Processes, CompiledProcesses import numpy as np, qsdsan.processes as pc, qsdsan as qs from qsdsan.utils import ospath, data_path +from qsdsan.processes._adm1 import ( + R, + create_adm1_cmps, + ADM1, + mass2mol_conversion, + T_correction_factor, + substr_inhibit, + non_compet_inhibit, + Hill_inhibit + ) # from scipy.optimize import brenth # from warnings import warn @@ -24,7 +34,6 @@ _path = ospath.join(data_path, 'process_data/_madm1.tsv') - #%% components # C_mw = get_mw({'C':1}) N_mw = get_mw({'N':1}) @@ -37,7 +46,7 @@ def create_madm1_cmps(set_thermo=True): # Components from the original ADM1 # ********************************* - _cmps = pc.create_adm1_cmps(False) + _cmps = create_adm1_cmps(False) S_aa = _cmps.S_aa X_pr = _cmps.X_pr S_aa.i_C = X_pr.i_C = 0.36890 @@ -161,10 +170,138 @@ def create_madm1_cmps(set_thermo=True): return cmps_madm1 #%% rate functions +{'S_su': 0, + 'S_aa': 1, + 'S_fa': 2, + 'S_va': 3, + 'S_bu': 4, + 'S_pro': 5, + 'S_ac': 6, + 'S_h2': 7, + 'S_ch4': 8, + 'S_IC': 9, + 'S_IN': 10, + 'S_IP': 11, + 'S_I': 12, + 'X_ch': 13, + 'X_pr': 14, + 'X_li': 15, + 'X_su': 16, + 'X_aa': 17, + 'X_fa': 18, + 'X_c4': 19, + 'X_pro': 20, + 'X_ac': 21, + 'X_h2': 22, + 'X_I': 23, + 'X_PHA': 24, + 'X_PP': 25, + 'X_PAO': 26, + 'S_K': 27, + 'S_Mg': 28, + 'S_SO4': 29, + 'S_IS': 30, + 'X_hSRB': 31, + 'X_aSRB': 32, + 'X_pSRB': 33, + 'X_c4SRB': 34, + 'S_S0': 35, + 'S_Fe3': 36, + 'S_Fe2': 37, + 'H2O': 38} + +def calc_pH(): + pass -def rhos_madm1(): +def calc_biogas(): pass +rhos = np.zeros(40) # 36 biochemical processes + 4 gas transfer processes +Cs = np.empty(36) + +def rhos_madm1(state_arr, params): + ks = params['rate_constants'] + Ks = params['half_sat_coeffs'] + K_PP = params['K_PP'] + K_so4 = params['K_so4'] + cmps = params['components'] + # n = len(cmps) + pH_LLs, pH_ULs = params['pH_limits'] + KS_IN = params['KS_IN'] + KS_IP = params['KS_IP'] + KI_nh3 = params['KI_nh3'] + KIs_h2 = params['KIs_h2'] + KIs_h2s = params['KIs_h2s'] + KHb = params['K_H_base'] + Kab = params['Ka_base'] + KH_dH = params['K_H_dH'] + Ka_dH = params['Ka_dH'] + kLa = params['kLa'] + T_base = params['T_base'] + + Cs[:7] = state_arr[13:20] # original ADM1 processes + Cs[7:11] = state_arr[19:23] + Cs[11:18] = state_arr[16:23] + Cs[18:23] = X_PAO = state_arr[26] # P extension processes + Cs[23:25] = X_PP, X_PHA = state_arr[[25,24]] + Cs[25:27] = state_arr[31] # S extension processes + Cs[27:29] = state_arr[32] + Cs[29:31] = state_arr[33] + Cs[31:34] = state_arr[34] + Cs[34:36] = state_arr[36] # Fe extension processes + + rhos[:-4] = ks * Cs + primary_substrates = state_arr[:8] + + rhos[3:11] *= substr_inhibit(primary_substrates, Ks[:8]) + c4 = primary_substrates[[3,4]] + if sum(c4) > 0: rhos[[6,7]] *= c4/sum(c4) + + vfas = primary_substrates[3:7] + rhos[18:22] *= substr_inhibit(vfas, Ks[8]) + if sum(vfas) > 0: rhos[18:22] *= vfas/sum(vfas) + if X_PAO > 0: rhos[18:22] *= substr_inhibit(X_PP/X_PAO, K_PP) + + srb_subs = np.flip(primary_substrates[3:]) + S_SO4, S_IS = state_arr[29:31] + rhos[[25,27,29,31,32]] *= substr_inhibit(srb_subs, Ks[-4:]) * substr_inhibit(S_SO4, K_so4) + if sum(srb_subs[-2:]) > 0: rhos[[31,32]] *= srb_subs[-2:]/sum(srb_subs[-2:]) + + #!!! why divide by 16 or 64? + S_h2 = primary_substrates[-1] + rhos[34] *= S_h2 / 16 + rhos[35] *= S_IS / 64 + +# ============================================================================= +# inhibition factors +# ============================================================================= + + S_IN, S_IP = state_arr[[10,11]] + I_nutrients = substr_inhibit(S_IN, KS_IN) * substr_inhibit(S_IP, KS_IP) + rhos[3:11] *= I_nutrients + rhos[[25,27,29,31,32]] *= I_nutrients + +# ============================================================================= +# #!!! place holder for pH +# ============================================================================= + pH, nh3 = calc_pH(state_arr, params) + Is_pH = Hill_inhibit(10**(-pH), pH_ULs, pH_LLs) + rhos[3:9] *= Is_pH[0] + rhos[9:11] *= Is_pH[1:3] + rhos[[25,27]] *= Is_pH[3:5] + rhos[[29,31,32]] *= Is_pH[-1] + + Is_h2 = non_compet_inhibit(S_h2, KIs_h2) + rhos[5:9] *= Is_h2 + Inh3 = non_compet_inhibit(nh3, KI_nh3) + rhos[9] *= Inh3 + + Z_h2s = calc_biogas() + Is_h2s = non_compet_inhibit(Z_h2s, KIs_h2s) + rhos[6:11] *= Is_h2s[:5] + rhos[[25,27,29,31,32]] *= Is_h2s[5:] + + #%% modified ADM1 class _load_components = settings.get_default_chemicals @@ -323,7 +460,7 @@ class ModifiedADM1(CompiledProcesses): _cmp_dependent_stoichio = ('K_XPP', 'Mg_XPP', 'MW_S0', 'MW_IS', 'i_mass_S0', 'i_mass_IS', 'i_mass_Fe2') - _stoichio_params = (*pc.ADM1._stoichio_params[5:], + _stoichio_params = (*ADM1._stoichio_params[5:], 'f_ch_xb', 'f_pr_xb', 'f_li_xb', 'f_xI_xb', 'f_sI_xb', 'f_va_pha', 'f_bu_pha', 'f_pro_pha', 'f_ac_pha', 'f_is_pro', 'f_is_bu', 'f_is_va', @@ -331,13 +468,14 @@ class ModifiedADM1(CompiledProcesses): *_cmp_dependent_stoichio ) _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'K_PP', 'K_so4', - 'pH_limits', 'KS_IN', 'KI_nh3', 'KIs_h2', + 'pH_limits', 'KS_IN', 'KS_IP', 'KI_nh3', 'KIs_h2', 'KIs_h2s' 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', - 'T_base', 'components', 'root', + 'T_base', 'components', + # 'root' ) - _acid_base_pairs = pc.ADM1._acid_base_pairs - _biogas_IDs = (*pc.ADM1._biogas_IDs, 'S_IS') - _biomass_IDs = (*pc.ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') + _acid_base_pairs = ADM1._acid_base_pairs + _biogas_IDs = (*ADM1._biogas_IDs, 'S_IS') + _biomass_IDs = (*ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') _T_base = 298.15 _K_H_base = [7.8e-4, 1.4e-3, 3.5e-2, 0.105] # biogas species Henry's Law constant [M/bar] _K_H_dH = [-4180, -14240, -19410, -19180] # Heat of reaction of liquid-gas transfer of biogas species [J/mol] From 4847f2f8ddc8908bb2ac2c522a3170803264e267 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 10 Oct 2023 15:57:26 -0700 Subject: [PATCH 07/16] MMP components and stoichiometry --- qsdsan/data/process_data/_madm1.tsv | 84 +++++++------ qsdsan/processes/_madm1.py | 189 +++++++++++++++++++++------- 2 files changed, 188 insertions(+), 85 deletions(-) diff --git a/qsdsan/data/process_data/_madm1.tsv b/qsdsan/data/process_data/_madm1.tsv index 2e124782..7493add1 100644 --- a/qsdsan/data/process_data/_madm1.tsv +++ b/qsdsan/data/process_data/_madm1.tsv @@ -1,37 +1,47 @@ - S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg S_SO4 S_IS X_hSRB X_aSRB X_pSRB X_c4SRB S_S0 S_Fe3 S_Fe2 -hydrolysis_carbs 1 ? ? ? -1 -hydrolysis_proteins 1 ? ? ? -1 -hydrolysis_lipids 1-f_fa_li f_fa_li ? ? ? -1 -uptake_sugars -1 (1-Y_su)*f_bu_su (1-Y_su)*f_pro_su (1-Y_su)*f_ac_su (1-Y_su)*f_h2_su ? ? ? Y_su -uptake_amino_acids -1 (1-Y_aa)*f_va_aa (1-Y_aa)*f_bu_aa (1-Y_aa)*f_pro_aa (1-Y_aa)*f_ac_aa (1-Y_aa)*f_h2_aa ? ? ? Y_aa -uptake_LCFA -1 (1-Y_fa)*f_ac_fa (1-Y_fa)*f_h2_fa ? ? ? Y_fa -uptake_valerate -1 (1-Y_c4)*f_pro_va (1-Y_c4)*f_ac_va (1-Y_c4)*f_h2_va ? ? ? Y_c4 -uptake_butyrate -1 (1-Y_c4)*f_ac_bu (1-Y_c4)*f_h2_bu ? ? ? Y_c4 -uptake_propionate -1 (1-Y_pro)*f_ac_pro (1-Y_pro)*f_h2_pro ? ? ? Y_pro -uptake_acetate -1 1-Y_ac ? ? ? Y_ac -uptake_h2 -1 1-Y_h2 ? ? ? Y_h2 -decay_Xsu ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xaa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xfa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xc4 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xpro ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xac ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -decay_Xh2 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb -storage_PHA_valerate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP -storage_PHA_butyrate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP -storage_PHA_propionate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP -storage_PHA_actate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP -lysis_XPAO ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 -lysis_XPP ? ? ? -1 K_XPP Mg_XPP -lysis_XPHA f_va_pha f_bu_pha f_pro_pha f_ac_pha ? ? ? -1 -growth_SRB_h2 -1 ? ? ? -1*(1-Y_hSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_hSRB Y_hSRB -decay_XhSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 -growth_SRB_acetate -1 ? ? ? -1*(1-Y_aSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_aSRB Y_aSRB -decay_XaSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 -growth_SRB_propionate -1 (1-Y_pSRB)*f_ac_pro ? ? ? -1*(1-Y_pSRB)*f_is_pro*i_mass_IS*(MW_S0/MW_IS) (1-Y_pSRB)*f_is_pro Y_pSRB -decay_XpSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 -growth_SRB_butyrate -1 (1-Y_c4SRB)*f_ac_bu ? ? ? -1*(1-Y_c4SRB)*f_is_bu*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_bu Y_c4SRB -growth_SRB_valerate -1 (1-Y_c4SRB)*f_pro_va (1-Y_c4SRB)*f_ac_va ? ? ? -1*(1-Y_c4SRB)*f_is_va*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_va Y_c4SRB -decay_Xc4SRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 -Fe3_reduction_h2 -1 -1*i_mass_Fe2 1 -Fe3_reduction_IS -1 (MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1) 1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) + S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg S_SO4 S_IS X_hSRB X_aSRB X_pSRB X_c4SRB S_S0 S_Fe3 S_Fe2 X_HFO_H X_HFO_L X_HFO_old X_HFO_HP X_HFO_LP X_HFO_HP_old X_HFO_LP_old +hydrolysis_carbs 1 ? ? ? -1 +hydrolysis_proteins 1 ? ? ? -1 +hydrolysis_lipids 1-f_fa_li f_fa_li ? ? ? -1 +uptake_sugars -1 (1-Y_su)*f_bu_su (1-Y_su)*f_pro_su (1-Y_su)*f_ac_su (1-Y_su)*f_h2_su ? ? ? Y_su +uptake_amino_acids -1 (1-Y_aa)*f_va_aa (1-Y_aa)*f_bu_aa (1-Y_aa)*f_pro_aa (1-Y_aa)*f_ac_aa (1-Y_aa)*f_h2_aa ? ? ? Y_aa +uptake_LCFA -1 (1-Y_fa)*f_ac_fa (1-Y_fa)*f_h2_fa ? ? ? Y_fa +uptake_valerate -1 (1-Y_c4)*f_pro_va (1-Y_c4)*f_ac_va (1-Y_c4)*f_h2_va ? ? ? Y_c4 +uptake_butyrate -1 (1-Y_c4)*f_ac_bu (1-Y_c4)*f_h2_bu ? ? ? Y_c4 +uptake_propionate -1 (1-Y_pro)*f_ac_pro (1-Y_pro)*f_h2_pro ? ? ? Y_pro +uptake_acetate -1 1-Y_ac ? ? ? Y_ac +uptake_h2 -1 1-Y_h2 ? ? ? Y_h2 +decay_Xsu ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xaa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xfa ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xc4 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xpro ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xac ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +decay_Xh2 ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb -1 f_xI_xb +storage_PHA_valerate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_butyrate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_propionate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +storage_PHA_actate -1 ? ? ? 1 -Y_PO4 Y_PO4*K_XPP Y_PO4*Mg_XPP +lysis_XPAO ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +lysis_XPP ? ? ? -1 K_XPP Mg_XPP +lysis_XPHA f_va_pha f_bu_pha f_pro_pha f_ac_pha ? ? ? -1 +growth_SRB_h2 -1 ? ? ? -1*(1-Y_hSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_hSRB Y_hSRB +decay_XhSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_acetate -1 ? ? ? -1*(1-Y_aSRB)*i_mass_IS*(MW_S0/MW_IS) 1-Y_aSRB Y_aSRB +decay_XaSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_propionate -1 (1-Y_pSRB)*f_ac_pro ? ? ? -1*(1-Y_pSRB)*f_is_pro*i_mass_IS*(MW_S0/MW_IS) (1-Y_pSRB)*f_is_pro Y_pSRB +decay_XpSRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +growth_SRB_butyrate -1 (1-Y_c4SRB)*f_ac_bu ? ? ? -1*(1-Y_c4SRB)*f_is_bu*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_bu Y_c4SRB +growth_SRB_valerate -1 (1-Y_c4SRB)*f_pro_va (1-Y_c4SRB)*f_ac_va ? ? ? -1*(1-Y_c4SRB)*f_is_va*i_mass_IS*(MW_S0/MW_IS) (1-Y_c4SRB)*f_is_va Y_c4SRB +decay_Xc4SRB ? ? ? f_sI_xb f_ch_xb f_pr_xb f_li_xb f_xI_xb -1 +reduction_HFO_H_h2 -1 1 -1*i_mass_Fe2 +reduction_HFO_L_h2 -1 1 -1*i_mass_Fe2 +reduction_HFO_H_IS -1 (MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) 1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1) +reduction_HFO_L_IS -1 (MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) 1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS) i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1) +aging_HFO_H -1 1 +aging_HFO_L -1 1 +fast_P_binding ? -1 1 +slow_P_sorption ? -1 1 +aging_HFO_HP -1 1 +aging_HFO_LP -1 1 +dissolution_HFO_HP ? 1 -1 +dissolution_HFO_LP ? 1 -1 diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index b8f83db6..537f0d6f 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -42,7 +42,7 @@ Fe_mw = get_mw({'Fe':1}) O_mw = get_mw({'O':1}) -def create_madm1_cmps(set_thermo=True): +def create_madm1_cmps(set_thermo=True, ASF_L=0.31, ASF_H=1.2): # Components from the original ADM1 # ********************************* @@ -101,24 +101,17 @@ def create_madm1_cmps(set_thermo=True): # Additional components for P, S, Fe extensions # ********************************************* S_IP = asm_cmps.S_PO4.copy('S_IP') - S_K = Component.from_chemical('S_K', chemical='K+', - description='Potassium ion', - measured_as='K', - particle_size='Soluble', - degradability='Undegradable', - organic=False) - S_Mg = Component.from_chemical('S_Mg', chemical='Mg2+', - description='Magnesium ion', - measured_as='Mg', - particle_size='Soluble', - degradability='Undegradable', - organic=False) - S_SO4 = Component.from_chemical('S_SO4', chemical='SO4-2', - description='Sulfate', - measured_as='S', - particle_size='Soluble', - degradability='Undegradable', - organic=False) + + ion_properties = dict( + particle_size='Soluble', + degradability='Undegradable', + organic=False) + S_K = Component.from_chemical('S_K', chemical='K+', description='Potassium ion', + measured_as='K', **ion_properties) + S_Mg = Component.from_chemical('S_Mg', chemical='Mg2+', description='Magnesium ion', + measured_as='Mg',**ion_properties) + S_SO4 = Component.from_chemical('S_SO4', chemical='SO4-2', description='Sulfate', + measured_as='S', **ion_properties) S_IS = Component.from_chemical('S_IS', chemical='H2S', description='Hydrogen sulfide', measured_as='COD', @@ -136,34 +129,81 @@ def create_madm1_cmps(set_thermo=True): X_c4SRB.description = 'sulfate-reducing biomass, utilizing butyrate and valerate' S_S0 = Component.from_chemical('S_S0', chemical='S', - description='Sulfur', + description='Elemental sulfur', measured_as='COD', particle_size='Soluble', degradability='Undegradable', organic=False) - S_Fe3 = Component.from_chemical('S_Fe3', chemical='Fe3+', - description='Iron (III)', - measured_as='Fe', - particle_size='Soluble', - degradability='Undegradable', - organic=False) - S_Fe2 = Component.from_chemical('S_Fe2', chemical='Fe2+', - description='Iron (II)', - measured_as='Fe', - particle_size='Soluble', - degradability='Undegradable', - organic=False) + S_Fe3 = Component.from_chemical('S_Fe3', chemical='Fe3+', description='Iron (III)', + measured_as='Fe',**ion_properties) + S_Fe2 = Component.from_chemical('S_Fe2', chemical='Fe2+', description='Iron (II)', + measured_as='Fe',**ion_properties) S_Fe2.i_COD = 0.5*O_mw/Fe_mw S_Fe2.measured_as = 'COD' + # Multiple mineral precipitation + # ****************************** + mineral_properties = dict( + particle_size='Particulate', + degradability='Undegradable', + organic=False) + + X_HFO_H = Component('X_HFO_H', formula='FeO(OH)', + description='Hydrous ferric oxide with high number of active sites', + measured_as='Fe',**mineral_properties) + X_HFO_L = X_HFO_H.copy('X_HFO_L') + X_HFO_L.description = 'Hydrous ferric oxide with low number of active sites' + + X_HFO_old = X_HFO_H.copy('X_HFO_old') + X_HFO_old.description = 'Inactive hydrous ferric oxide' + + X_HFO_HP = Component('X_HFO_HP', formula=f'FeO(OH)P{ASF_H}', + description='X_HFO_H with phosphorus-bounded adsorption sites', + measured_as='Fe', **mineral_properties) + X_HFO_HP_old = X_HFO_HP.copy('X_HFO_HP_old') + X_HFO_HP_old.description = 'Old ' + X_HFO_HP.description + + X_HFO_LP = Component('X_HFO_LP', formula=f'FeO(OH)P{ASF_L}', + description='X_HFO_L with phosphorus-bounded adsorption sites', + measured_as='Fe', **mineral_properties) + X_HFO_LP_old = X_HFO_LP.copy('X_HFO_LP_old') + X_HFO_LP_old.description = 'Old ' + X_HFO_LP.description + + X_CCM = Component.from_chemical('X_CCM', chemical='calcite', description='Calcite', **mineral_properties) + X_ACC = Component.from_chemical('X_ACC', chemical='aragonite', description='Aragonite', **mineral_properties) + X_ACP = Component.from_chemical('X_ACP', chemical='Ca3(PO4)2', description='Amorphous calcium phosphate', **mineral_properties) + X_HAP = Component.from_chemical('X_HAP', chemical='hydroxylapatite', description='Hydroxylapatite', **mineral_properties) + X_DCPD = Component.from_chemical('X_DCPD', chemical='CaHPO4', description='Dicalcium phosphate', **mineral_properties) + X_OCP = Component('X_OCP', formula='Ca4HP3O12', description='Octacalcium phosphate', **mineral_properties) + X_struv = Component.from_chemical('X_struv', chemical='MgNH4PO4', description='Struvite', **mineral_properties) + X_newb = Component.from_chemical('X_newb', chemical='MgHPO4', description='Newberyite', **mineral_properties) + X_magn = Component.from_chemical('X_magn', chemical='MgCO3', description='Magnesite', **mineral_properties) + X_kstruv = Component('X_kstruv', formula='MgKPO4', description='K-struvite', **mineral_properties) + X_FeS = Component.from_chemical('X_FeS', chemical='FeS', description='Iron sulfide', **mineral_properties) + X_Fe3PO42 = Component('X_Fe3PO42', formula='Fe3(PO4)2', description='Ferrous phosphate', **mineral_properties) + X_AlPO4 = Component.from_chemical('X_AlPO4', chemical='AlPO4', description='Aluminum phosphate', **mineral_properties) + + S_Ca = Component.from_chemical('S_Ca', chemical='Ca2+', description='Calsium ion', + measured_as='Ca', **ion_properties) + S_Al = Component.from_chemical('S_Al', chemical='Al3+', description='Aluminum ion', + measured_as='Al', **ion_properties) + S_Na = Component.from_chemical('S_Na', chemical='Na+', description='Sodium ion', + measured_as='Na', **ion_properties) + S_Cl = Component.from_chemical('S_Cl', chemical='Cl-', description='Chloride', + measured_as='Cl', **ion_properties) + cmps_madm1 = Components([_cmps.S_su, S_aa, S_fa, _cmps.S_va, S_bu, S_pro, S_ac, _cmps.S_h2, _cmps.S_ch4, _cmps.S_IC, _cmps.S_IN, S_IP, S_I, X_ch, X_pr, X_li, *adm1_biomass, X_I, X_PHA, asm_cmps.X_PP, X_PAO, S_K, S_Mg, S_SO4, S_IS, X_hSRB, X_aSRB, X_pSRB, X_c4SRB, - S_S0, S_Fe3, S_Fe2, - _cmps.H2O]) + S_S0, S_Fe3, S_Fe2, X_HFO_H, X_HFO_L, X_HFO_old, + X_HFO_HP, X_HFO_LP, X_HFO_HP_old, X_HFO_LP_old, + S_Ca, S_Al, X_CCM, X_ACC, X_ACP, X_HAP, X_DCPD, + X_OCP, X_struv, X_newb, X_magn, X_kstruv, X_FeS, + X_Fe3PO42, X_AlPO4, + S_Na, S_Cl, _cmps.H2O]) cmps_madm1.default_compile() if set_thermo: qs.set_thermo(cmps_madm1) @@ -296,7 +336,10 @@ def rhos_madm1(state_arr, params): Inh3 = non_compet_inhibit(nh3, KI_nh3) rhos[9] *= Inh3 - Z_h2s = calc_biogas() +# ============================================================================= +# !!! place holder for gas-liquid transfer +# ============================================================================= + Z_h2s = calc_biogas() # should be a function of pH Is_h2s = non_compet_inhibit(Z_h2s, KIs_h2s) rhos[6:11] *= Is_h2s[:5] rhos[[25,27,29,31,32]] *= Is_h2s[5:] @@ -309,7 +352,7 @@ def rhos_madm1(state_arr, params): @chemicals_user class ModifiedADM1(CompiledProcesses): """ - Modified Anaerobic Digestion Model no.1 [1]_ + Modified Anaerobic Digestion Model no.1 [1]_, [2]_ Parameters ---------- @@ -446,11 +489,17 @@ class ModifiedADM1(CompiledProcesses): References ---------- .. [1] Flores-Alsina, X., Solon, K., Kazadi Mbamba, C., Tait, S., - Gernaey, K. V., Jeppsson, U., ; Batstone, D. J. (2016). + Gernaey, K. V., Jeppsson, U., Batstone, D. J. (2016). Modelling phosphorus (P), sulfur (S) and iron (Fe) interactions for dynamic simulations of anaerobic digestion processes. Water Research, 95, 370–382. https://doi.org/10.1016/J.WATRES.2016.03.012 - + .. [2] Solon, K., Flores-Alsina, X., Kazadi Mbamba, C., Ikumi, D., + Volcke, E. I. P., Vaneeckhaute, C., Ekama, G., Vanrolleghem, P. A., + Batstone, D. J., Gernaey, K. v., Jeppsson, U. (2017). Plant-wide + modelling of phosphorus transformations in wastewater treatment systems: + Impacts of control and operational strategies. Water Research, + 113, 97–110. https://doi.org/10.1016/J.WATRES.2017.02.007 + See Also -------- `qsdsan.processes.ADM1 `_ @@ -476,6 +525,9 @@ class ModifiedADM1(CompiledProcesses): _acid_base_pairs = ADM1._acid_base_pairs _biogas_IDs = (*ADM1._biogas_IDs, 'S_IS') _biomass_IDs = (*ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') + _precipitate_IDs = ('X_CCM', 'X_ACC', 'X_ACP', 'X_HAP', 'X_DCPD', 'X_OCP', + 'X_struv', 'X_newb', 'X_magn', 'X_kstruv', + 'X_FeS', 'X_Fe3PO42', 'X_AlPO4') _T_base = 298.15 _K_H_base = [7.8e-4, 1.4e-3, 3.5e-2, 0.105] # biogas species Henry's Law constant [M/bar] _K_H_dH = [-4180, -14240, -19410, -19180] # Heat of reaction of liquid-gas transfer of biogas species [J/mol] @@ -515,6 +567,47 @@ def __new__(cls, components=None, path=None, parameters=cls._stoichio_params, compile=False) + precipitation = [] + for i in cls._precipitate_IDs[:-3]: + new_p = Process('precipitation_%s' % i.lstrip('X_'), + reaction='[?]S_IC + [?]S_IN + [?]S_IP + [?]S_K + [?]S_Mg + [?]S_Ca -> %s' % i, + ref_component=i, + conserved_for=('C', 'N', 'P', 'K', 'Mg', 'Ca'), + parameters=()) + precipitation.append(new_p) + + i_mass_IS = cmps.S_IS.i_mass + i_mass_Fe2 = cmps.S_Fe2.i_mass + FeS_mw = cmps.X_FeS.chem_MW + new_p = Process('precipitation_FeS', + reaction={'S_Fe2': -Fe_mw/FeS_mw/i_mass_Fe2, + 'S_IS': -S_mw/FeS_mw/i_mass_IS, + 'X_FeS': 1}, + ref_component='X_FeS', + conserved_for=()) + precipitation.append(new_p) + + Fe3PO42_mw = cmps.X_Fe3PO42.chem_MW + new_p = Process('precipitation_Fe3PO42', + reaction={'S_Fe2': -3*Fe_mw/Fe3PO42_mw/i_mass_Fe2, + 'S_IP': '?', + 'X_Fe3PO42': 1}, + ref_component='X_Fe3PO42', + conserved_for=('P',)) + precipitation.append(new_p) + + AlPO4_mw = cmps.X_AlPO4.chem_MW + Al_mw = cmps.S_Al.chem_MW + new_p = Process('precipitation_AlPO4', + reaction={'S_Al': -Al_mw/AlPO4_mw, + 'S_IP': '?', + 'X_AlPO4': 1}, + ref_component='X_AlPO4', + conserved_for=('P',)) + precipitation.append(new_p) + + self.extend(precipitation) + gas_transfer = [] for i in cls._biogas_IDs: new_p = Process('%s_transfer' % i.lstrip('S_'), @@ -537,7 +630,7 @@ def __new__(cls, components=None, path=None, Y_PO4, Y_hSRB, Y_aSRB, Y_pSRB, Y_c4SRB, cmps.X_PP.i_K, cmps.X_PP.i_Mg, cmps.S_S0.chem_MW, cmps.S_IS.chem_MW, - cmps.S_S0.i_mass, cmps.S_IS.i_mass, cmps.S_Fe2.i_mass) + cmps.S_S0.i_mass, i_mass_IS, i_mass_Fe2) pH_limits = np.array([pH_limits_aa, pH_limits_ac, pH_limits_h2, pH_limits_h2_SRB, pH_limits_ac_SRB, pH_limits_aa_SRB]).T @@ -565,14 +658,14 @@ def __new__(cls, components=None, path=None, dct = self.__dict__ dct.update(kwargs) - self.set_rate_function(rhos_madm1) + # self.set_rate_function(rhos_madm1) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) - self.rate_function._params = dict(zip(cls._kinetic_params, - [ks, Ks, K_PP, K_so4, - pH_limits, KS_IN*N_mw, KS_IP*P_mw, - KI_nh3, KIs_h2, KIs_h2s, - Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, - cls.T_base, self._components, - # root, - ])) + # self.rate_function._params = dict(zip(cls._kinetic_params, + # [ks, Ks, K_PP, K_so4, + # pH_limits, KS_IN*N_mw, KS_IP*P_mw, + # KI_nh3, KIs_h2, KIs_h2s, + # Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, + # cls.T_base, self._components, + # # root, + # ])) return self \ No newline at end of file From ff6045ae51aa84addffed746f6903e6576f22e9a Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 10 Oct 2023 16:24:47 -0700 Subject: [PATCH 08/16] additional references --- qsdsan/processes/_madm1.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 537f0d6f..cd34e685 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -352,7 +352,7 @@ def rhos_madm1(state_arr, params): @chemicals_user class ModifiedADM1(CompiledProcesses): """ - Modified Anaerobic Digestion Model no.1 [1]_, [2]_ + Modified Anaerobic Digestion Model no.1 [1]_, [2]_, [3]_ Parameters ---------- @@ -498,7 +498,10 @@ class ModifiedADM1(CompiledProcesses): Batstone, D. J., Gernaey, K. v., Jeppsson, U. (2017). Plant-wide modelling of phosphorus transformations in wastewater treatment systems: Impacts of control and operational strategies. Water Research, - 113, 97–110. https://doi.org/10.1016/J.WATRES.2017.02.007 + 113, 97–110. https://doi.org/10.1016/J.WATRES.2017.02.007 + .. [3] Hauduc, H., Takács, I., Smith, S., Szabo, A., Murthy, S., Daigger, G. T., + Spérandio, M. (2015). A dynamic physicochemical model for chemical phosphorus + removal. Water Research, 73, 157–170. https://doi.org/10.1016/J.WATRES.2014.12.053 See Also -------- @@ -658,8 +661,8 @@ def __new__(cls, components=None, path=None, dct = self.__dict__ dct.update(kwargs) - # self.set_rate_function(rhos_madm1) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) + # self.set_rate_function(rhos_madm1) # self.rate_function._params = dict(zip(cls._kinetic_params, # [ks, Ks, K_PP, K_so4, # pH_limits, KS_IN*N_mw, KS_IP*P_mw, From 6731162ff13ff36a16be1e94006e6eaa7bc88fea Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 11 Oct 2023 11:13:24 -0700 Subject: [PATCH 09/16] updated Fe extension with HFO --- qsdsan/processes/_madm1.py | 107 ++++++++++++++++++------------------- 1 file changed, 52 insertions(+), 55 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index cd34e685..66038651 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -43,6 +43,23 @@ O_mw = get_mw({'O':1}) def create_madm1_cmps(set_thermo=True, ASF_L=0.31, ASF_H=1.2): + ''' + Create a set of components for the modified ADM1. + + Parameters + ---------- + set_thermo : bool, optional + Whether to set thermo with the returned set of components. The default is True. + ASF_L : float, optional + Active site factor for X_HFO_L [mol P sites/mol Fe]. The default is 0.31. + ASF_H : float, optional + Active site factor for X_HFO_H [mol P sites/mol Fe]. The default is 1.2. + + Returns + ------- + cmps_madm1 : class:`CompiledComponents` + + ''' # Components from the original ADM1 # ********************************* @@ -210,45 +227,23 @@ def create_madm1_cmps(set_thermo=True, ASF_L=0.31, ASF_H=1.2): return cmps_madm1 #%% rate functions -{'S_su': 0, - 'S_aa': 1, - 'S_fa': 2, - 'S_va': 3, - 'S_bu': 4, - 'S_pro': 5, - 'S_ac': 6, - 'S_h2': 7, - 'S_ch4': 8, - 'S_IC': 9, - 'S_IN': 10, - 'S_IP': 11, - 'S_I': 12, - 'X_ch': 13, - 'X_pr': 14, - 'X_li': 15, - 'X_su': 16, - 'X_aa': 17, - 'X_fa': 18, - 'X_c4': 19, - 'X_pro': 20, - 'X_ac': 21, - 'X_h2': 22, - 'X_I': 23, - 'X_PHA': 24, - 'X_PP': 25, - 'X_PAO': 26, - 'S_K': 27, - 'S_Mg': 28, - 'S_SO4': 29, - 'S_IS': 30, - 'X_hSRB': 31, - 'X_aSRB': 32, - 'X_pSRB': 33, - 'X_c4SRB': 34, - 'S_S0': 35, - 'S_Fe3': 36, - 'S_Fe2': 37, - 'H2O': 38} + +"https://wiki.dynamita.com/en/biokinetic_process_models#chemical-phosphorus-removal-with-metal-salts-addition-iron-or-aluminium" + +{'S_su': 0, 'S_aa': 1, 'S_fa': 2, 'S_va': 3, 'S_bu': 4, 'S_pro': 5, 'S_ac': 6, 'S_h2': 7, + 'S_ch4': 8, 'S_IC': 9, 'S_IN': 10, 'S_IP': 11, 'S_I': 12, + 'X_ch': 13, 'X_pr': 14, 'X_li': 15, + 'X_su': 16, 'X_aa': 17, 'X_fa': 18, 'X_c4': 19, 'X_pro': 20, 'X_ac': 21, 'X_h2': 22, 'X_I': 23, + 'X_PHA': 24, 'X_PP': 25, 'X_PAO': 26, 'S_K': 27, 'S_Mg': 28, + 'S_SO4': 29, 'S_IS': 30, 'X_hSRB': 31, 'X_aSRB': 32, 'X_pSRB': 33, 'X_c4SRB': 34, + 'S_S0': 35, 'S_Fe3': 36, 'S_Fe2': 37, + 'X_HFO_H': 38, 'X_HFO_L': 39, 'X_HFO_old': 40, 'X_HFO_HP': 41, 'X_HFO_LP': 42, 'X_HFO_HP_old': 43, 'X_HFO_LP_old': 44, + 'S_Ca': 45, 'S_Al': 46, + 'X_CCM': 47, 'X_ACC': 48, 'X_ACP': 49, 'X_HAP': 50, 'X_DCPD': 51, 'X_OCP': 52, + 'X_struv': 53, 'X_newb': 54, 'X_magn': 55, 'X_kstruv': 56, + 'X_FeS': 57, 'X_Fe3PO42': 58, + 'X_AlPO4': 59, + 'S_Na': 60, 'S_Cl': 61, 'H2O': 62} def calc_pH(): pass @@ -256,8 +251,8 @@ def calc_pH(): def calc_biogas(): pass -rhos = np.zeros(40) # 36 biochemical processes + 4 gas transfer processes -Cs = np.empty(36) +rhos = np.zeros(38+8+13+4) # 38 biological + 8 chemical P removal by HFO + 13 MMP + 4 gas transfer +Cs = np.empty(38) def rhos_madm1(state_arr, params): ks = params['rate_constants'] @@ -288,9 +283,9 @@ def rhos_madm1(state_arr, params): Cs[27:29] = state_arr[32] Cs[29:31] = state_arr[33] Cs[31:34] = state_arr[34] - Cs[34:36] = state_arr[36] # Fe extension processes + Cs[34:36] = Cs[36:38] = state_arr[38:40] # Fe extension processes - rhos[:-4] = ks * Cs + rhos[:38] = ks * Cs primary_substrates = state_arr[:8] rhos[3:11] *= substr_inhibit(primary_substrates, Ks[:8]) @@ -309,12 +304,11 @@ def rhos_madm1(state_arr, params): #!!! why divide by 16 or 64? S_h2 = primary_substrates[-1] - rhos[34] *= S_h2 / 16 - rhos[35] *= S_IS / 64 + rhos[34:36] *= S_h2 / 16 + rhos[36:38] *= S_IS / 64 -# ============================================================================= -# inhibition factors -# ============================================================================= + # inhibition factors + # ****************** S_IN, S_IP = state_arr[[10,11]] I_nutrients = substr_inhibit(S_IN, KS_IN) * substr_inhibit(S_IP, KS_IP) @@ -322,7 +316,7 @@ def rhos_madm1(state_arr, params): rhos[[25,27,29,31,32]] *= I_nutrients # ============================================================================= -# #!!! place holder for pH +# !!! place holder for pH # ============================================================================= pH, nh3 = calc_pH(state_arr, params) Is_pH = Hill_inhibit(10**(-pH), pH_ULs, pH_LLs) @@ -339,7 +333,7 @@ def rhos_madm1(state_arr, params): # ============================================================================= # !!! place holder for gas-liquid transfer # ============================================================================= - Z_h2s = calc_biogas() # should be a function of pH + Z_h2s = calc_biogas() # should be a function of pH, like co2 and nh3 Is_h2s = non_compet_inhibit(Z_h2s, KIs_h2s) rhos[6:11] *= Is_h2s[:5] rhos[[25,27,29,31,32]] *= Is_h2s[5:] @@ -442,9 +436,12 @@ class ModifiedADM1(CompiledProcesses): K_so4_c4SRB : float, optional Sulfate half saturation coefficient of SRB uptaking butyrate or valerate [kg S/m3]. The default is 6.413e-3. - k_Fe3t2 : float, optional - Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)]. - The default is 1.79e7. + k_Fe3t2_h2 : float, optional + Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)] using hydrogen + as electron donor. The default is 1.79e7. + k_Fe3t2_is : float, optional + Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)] using sulfide + as electron donor. The default is 1.79e7. KS_IP : float, optional Inorganic phosphorus (nutrient) inhibition coefficient for soluble substrate uptake [M]. The default is 2e-5. @@ -552,7 +549,7 @@ def __new__(cls, components=None, path=None, b_hSRB=0.02, b_aSRB=0.02, b_pSRB=0.02, b_c4SRB=0.02, K_hSRB=5.96e-6, K_aSRB=0.176, K_pSRB=0.088, K_c4SRB=0.1739, K_so4_hSRB=1.04e-4*S_mw, K_so4_aSRB=2e-4*S_mw, K_so4_pSRB=2e-4*S_mw, K_so4_c4SRB=2e-4*S_mw, - k_Fe3t2=1e9/Fe_mw, + k_Fe3t2_h2=1e9/Fe_mw, k_Fe3t2_is=1e9/Fe_mw, KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, KS_IP=2e-5, KI_h2s_c4=0.481, KI_h2s_pro=0.481, KI_h2s_ac=0.460, KI_h2s_h2=0.400, KI_h2s_c4SRB=0.520, KI_h2s_pSRB=0.520, KI_h2s_aSRB=0.499, KI_h2s_hSRB=0.499, @@ -643,7 +640,7 @@ def __new__(cls, components=None, path=None, b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2, # original ADM1 q_pha, q_pha, q_pha, q_pha, b_pao, b_pp, b_pha, # P extension k_hSRB, b_hSRB, k_aSRB, b_aSRB, k_pSRB, b_pSRB, k_c4SRB, k_c4SRB, b_c4SRB, # S extension - k_Fe3t2, k_Fe3t2)) # Fe extension + k_Fe3t2_h2, k_Fe3t2_h2, k_Fe3t2_is, k_Fe3t2_is)) # Fe extension + HFO Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2, # original ADM1 K_A, # P extension From 3eeb04520a40a7f0a5eee71ac015020fafdbd2cd Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 11 Oct 2023 11:35:36 -0700 Subject: [PATCH 10/16] scale stoichiometric coefficients to match default rate equation --- qsdsan/_process.py | 5 +++-- qsdsan/processes/_madm1.py | 4 ++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/qsdsan/_process.py b/qsdsan/_process.py index df1a812b..35ae8ff1 100644 --- a/qsdsan/_process.py +++ b/qsdsan/_process.py @@ -706,8 +706,9 @@ def _normalize_stoichiometry(self, new_ref): self._stoichiometry = [v/factor for v in stoich] def _normalize_rate_eq(self, new_ref): - factor = abs(self._stoichiometry[self._components.index(str(new_ref))]) - self._rate_equation *= factor + if self._rate_equation: + factor = abs(self._stoichiometry[self._components.index(str(new_ref))]) + self._rate_equation *= factor def show(self): info = f"Process: {self.ID}" diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 66038651..47d06dd9 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -567,6 +567,10 @@ def __new__(cls, components=None, path=None, parameters=cls._stoichio_params, compile=False) + for i in ('fast_P_binding', 'slow_P_sorption', 'dissolution_HFO_HP', 'dissolution_HFO_LP'): + p = getattr(self, i) + p.ref_component = 'S_IP' + precipitation = [] for i in cls._precipitate_IDs[:-3]: new_p = Process('precipitation_%s' % i.lstrip('X_'), From e8821b13631c58de73e4aa30daa3ef334af0ee06 Mon Sep 17 00:00:00 2001 From: Yalin Date: Wed, 11 Oct 2023 20:25:09 -0400 Subject: [PATCH 11/16] allowing additional properties to be stored with `WasteStream` --- qsdsan/_waste_stream.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/qsdsan/_waste_stream.py b/qsdsan/_waste_stream.py index db79c729..ceeb6296 100644 --- a/qsdsan/_waste_stream.py +++ b/qsdsan/_waste_stream.py @@ -40,12 +40,13 @@ _defined_composite_vars = ('COD', 'BOD5', 'BOD', 'uBOD', 'NOD', 'ThOD', 'cnBOD', 'C', 'N', 'P', 'K', 'Mg', 'Ca', 'solids', 'charge') -_common_composite_vars = ('_COD', '_BOD', '_uBOD', '_TC', '_TOC', '_TN', +_common_composite_vars = ('_COD', '_BOD', '_uBOD', '_ThOD', '_cnBOD', + '_TC', '_TOC', '_TN', '_TKN', '_TP', '_TK', '_TMg', '_TCa', - '_dry_mass', '_charge', '_ThOD', '_cnBOD') + '_dry_mass', '_charge',) _ws_specific_slots = (*_common_composite_vars, - '_pH', '_SAlk', '_ratios', + '_pH', '_SAlk', '_ratios', 'additional_properties', # '_stream_impact_item', (pls keep this here, might be useful in debugging) '_state', '_dstate', '_scope') @@ -249,7 +250,9 @@ class WasteStream(SanStream): TO BE IMPLEMENTED stream_impact_item : :class:`StreamImpactItem` The :class:`StreamImpactItem` this stream is linked to. - component_flows : kwargs + additional_properties : dict + Additional properties (e.g., turbidity). + component_flows : dict Component flow data. @@ -271,20 +274,28 @@ def __init__(self, ID='', flow=(), phase='l', T=298.15, P=101325., pH=7., SAlk=2.5, COD=None, BOD=None, uBOD=None, ThOD=None, cnBOD=None, TC=None, TOC=None, TN=None, TKN=None, TP=None, TK=None, - TMg=None, TCa=None, dry_mass=None, charge=None, ratios=None, - stream_impact_item=None, **component_flows): + TMg=None, TCa=None, + dry_mass=None, charge=None, ratios=None, + stream_impact_item=None, additional_properties={}, + **component_flows): SanStream.__init__(self=self, ID=ID, flow=flow, phase=phase, T=T, P=P, units=units, price=price, thermo=thermo, stream_impact_item=stream_impact_item, **component_flows) - self._init_ws(pH, SAlk, COD, BOD, uBOD, TC, TOC, TN, TKN, - TP, TK, TMg, TCa, ThOD, cnBOD, dry_mass, charge, ratios) + self._init_ws(pH=pH, SAlk=SAlk, COD=COD, BOD=BOD, uBOD=uBOD, + ThOD=ThOD, cnBOD=cnBOD, + TC=TC, TOC=TOC, TN=TN, TKN=TKN, + TP=TP, TK=TK, TMg=TMg, TCa=TCa, + dry_mass=dry_mass, charge=charge, ratios=ratios, + additional_properties=additional_properties, + ) def _init_ws(self, pH=7., SAlk=None, COD=None, BOD=None, uBOD=None, ThOD=None, cnBOD=None, TC=None, TOC=None, TN=None, TKN=None, TP=None, TK=None, TMg=None, TCa=None, - dry_mass=None, charge=None, ratios=None): + dry_mass=None, charge=None, ratios=None, + additional_properties={}): self._pH = pH self._SAlk = SAlk @@ -306,6 +317,7 @@ def _init_ws(self, pH=7., SAlk=None, COD=None, BOD=None, self._ratios = ratios self._state = None self._dstate = None + self.additional_properties = {} @staticmethod def from_stream(stream, ID='', **kwargs): @@ -2241,6 +2253,11 @@ def dry_mass(self): def density(self): '''[float] Density of the stream, in mg/L (kg/m3).''' return 0. + + @property + def additional_property(self): + '''[dict] Additional properties (e.g., turbidity).''' + return {} #!!! Keep this up-to-date with WasteStream # @property From 56c45f282816aa571998dbe7f41744705def7c9c Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 11 Oct 2023 18:34:06 -0700 Subject: [PATCH 12/16] kinetics for multiple mineral precipitation --- qsdsan/processes/_adm1.py | 1 + qsdsan/processes/_madm1.py | 165 ++++++++++++++++++++++++++++--------- 2 files changed, 127 insertions(+), 39 deletions(-) diff --git a/qsdsan/processes/_adm1.py b/qsdsan/processes/_adm1.py index 32d99819..cb876244 100644 --- a/qsdsan/processes/_adm1.py +++ b/qsdsan/processes/_adm1.py @@ -190,6 +190,7 @@ def mass2mol_conversion(cmps): # return np.exp(theta * (T2-T1)) def T_correction_factor(T1, T2, delta_H): + if T1 == T2: return 1 return np.exp(delta_H/(R*100) * (1/T1 - 1/T2)) # R converted to SI # def calc_Kas(pKas, T_base, T_op, theta): diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index 47d06dd9..ade77197 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -228,22 +228,26 @@ def create_madm1_cmps(set_thermo=True, ASF_L=0.31, ASF_H=1.2): #%% rate functions -"https://wiki.dynamita.com/en/biokinetic_process_models#chemical-phosphorus-removal-with-metal-salts-addition-iron-or-aluminium" - -{'S_su': 0, 'S_aa': 1, 'S_fa': 2, 'S_va': 3, 'S_bu': 4, 'S_pro': 5, 'S_ac': 6, 'S_h2': 7, - 'S_ch4': 8, 'S_IC': 9, 'S_IN': 10, 'S_IP': 11, 'S_I': 12, - 'X_ch': 13, 'X_pr': 14, 'X_li': 15, - 'X_su': 16, 'X_aa': 17, 'X_fa': 18, 'X_c4': 19, 'X_pro': 20, 'X_ac': 21, 'X_h2': 22, 'X_I': 23, - 'X_PHA': 24, 'X_PP': 25, 'X_PAO': 26, 'S_K': 27, 'S_Mg': 28, - 'S_SO4': 29, 'S_IS': 30, 'X_hSRB': 31, 'X_aSRB': 32, 'X_pSRB': 33, 'X_c4SRB': 34, - 'S_S0': 35, 'S_Fe3': 36, 'S_Fe2': 37, - 'X_HFO_H': 38, 'X_HFO_L': 39, 'X_HFO_old': 40, 'X_HFO_HP': 41, 'X_HFO_LP': 42, 'X_HFO_HP_old': 43, 'X_HFO_LP_old': 44, - 'S_Ca': 45, 'S_Al': 46, - 'X_CCM': 47, 'X_ACC': 48, 'X_ACP': 49, 'X_HAP': 50, 'X_DCPD': 51, 'X_OCP': 52, - 'X_struv': 53, 'X_newb': 54, 'X_magn': 55, 'X_kstruv': 56, - 'X_FeS': 57, 'X_Fe3PO42': 58, - 'X_AlPO4': 59, - 'S_Na': 60, 'S_Cl': 61, 'H2O': 62} +# https://wiki.dynamita.com/en/biokinetic_process_models#chemical-phosphorus-removal-with-metal-salts-addition-iron-or-aluminium + +# ============================================================================= +# state_variable_indices = { +# 'S_su': 0, 'S_aa': 1, 'S_fa': 2, 'S_va': 3, 'S_bu': 4, 'S_pro': 5, 'S_ac': 6, 'S_h2': 7, +# 'S_ch4': 8, 'S_IC': 9, 'S_IN': 10, 'S_IP': 11, 'S_I': 12, +# 'X_ch': 13, 'X_pr': 14, 'X_li': 15, +# 'X_su': 16, 'X_aa': 17, 'X_fa': 18, 'X_c4': 19, 'X_pro': 20, 'X_ac': 21, 'X_h2': 22, 'X_I': 23, +# 'X_PHA': 24, 'X_PP': 25, 'X_PAO': 26, 'S_K': 27, 'S_Mg': 28, +# 'S_SO4': 29, 'S_IS': 30, 'X_hSRB': 31, 'X_aSRB': 32, 'X_pSRB': 33, 'X_c4SRB': 34, +# 'S_S0': 35, 'S_Fe3': 36, 'S_Fe2': 37, +# 'X_HFO_H': 38, 'X_HFO_L': 39, 'X_HFO_old': 40, 'X_HFO_HP': 41, 'X_HFO_LP': 42, 'X_HFO_HP_old': 43, 'X_HFO_LP_old': 44, +# 'S_Ca': 45, 'S_Al': 46, +# 'X_CCM': 47, 'X_ACC': 48, 'X_ACP': 49, 'X_HAP': 50, 'X_DCPD': 51, 'X_OCP': 52, +# 'X_struv': 53, 'X_newb': 54, 'X_magn': 55, 'X_kstruv': 56, +# 'X_FeS': 57, 'X_Fe3PO42': 58, +# 'X_AlPO4': 59, +# 'S_Na': 60, 'S_Cl': 61, 'H2O': 62 +# } +# ============================================================================= def calc_pH(): pass @@ -251,8 +255,13 @@ def calc_pH(): def calc_biogas(): pass +def pcm(): + pass + + rhos = np.zeros(38+8+13+4) # 38 biological + 8 chemical P removal by HFO + 13 MMP + 4 gas transfer -Cs = np.empty(38) +Cs = np.empty(38+8) +sum_stoichios = np.array([2, 2, 5, 9, 3, 8, 3, 3, 2, 3, 2, 2]) def rhos_madm1(state_arr, params): ks = params['rate_constants'] @@ -272,7 +281,10 @@ def rhos_madm1(state_arr, params): KH_dH = params['K_H_dH'] Ka_dH = params['Ka_dH'] kLa = params['kLa'] + k_cryst = params['k_cryst'] + n_cryst = params['n_cryst'] T_base = params['T_base'] + # T_temp = params.pop('T_temp', T_base) Cs[:7] = state_arr[13:20] # original ADM1 processes Cs[7:11] = state_arr[19:23] @@ -283,9 +295,10 @@ def rhos_madm1(state_arr, params): Cs[27:29] = state_arr[32] Cs[29:31] = state_arr[33] Cs[31:34] = state_arr[34] - Cs[34:36] = Cs[36:38] = state_arr[38:40] # Fe extension processes + Cs[34:36] = Cs[36:38] = Cs[38:40] = Cs[40:42] = state_arr[38:40] # Fe extension processes + HFO module + Cs[42:44] = Cs[44:46] = state_arr[41:43] - rhos[:38] = ks * Cs + rhos[:46] = ks * Cs primary_substrates = state_arr[:8] rhos[3:11] *= substr_inhibit(primary_substrates, Ks[:8]) @@ -299,7 +312,7 @@ def rhos_madm1(state_arr, params): srb_subs = np.flip(primary_substrates[3:]) S_SO4, S_IS = state_arr[29:31] - rhos[[25,27,29,31,32]] *= substr_inhibit(srb_subs, Ks[-4:]) * substr_inhibit(S_SO4, K_so4) + rhos[[25,27,29,31,32]] *= substr_inhibit(srb_subs, Ks[9:13]) * substr_inhibit(S_SO4, K_so4) if sum(srb_subs[-2:]) > 0: rhos[[31,32]] *= srb_subs[-2:]/sum(srb_subs[-2:]) #!!! why divide by 16 or 64? @@ -307,6 +320,11 @@ def rhos_madm1(state_arr, params): rhos[34:36] *= S_h2 / 16 rhos[36:38] *= S_IS / 64 + KPbind, KPdiss = Ks[-2:] + S_IP = state_arr[11] + rhos[40:42] *= substr_inhibit(S_IP, KPbind) + rhos[44:46] *= non_compet_inhibit(S_IP, KPdiss) + # inhibition factors # ****************** @@ -338,11 +356,31 @@ def rhos_madm1(state_arr, params): rhos[6:11] *= Is_h2s[:5] rhos[[25,27,29,31,32]] *= Is_h2s[5:] - +# ============================================================================= +# !!! place holder for PCM (speciation) +# ============================================================================= + acts = pcm(state_arr) + SIs = np.maximum(1.0, saturation_index(acts)) # should be an array + rhos[46:59] = k_cryst * state_arr[47:60] * (SIs**(1/sum_stoichios) - 1)**n_cryst #%% modified ADM1 class _load_components = settings.get_default_chemicals +def fun(q_aging_H=450.0, q_aging_L=0.1, q_Pcoprec=360, q_Pbinding=0.3, q_diss_H=36.0, q_diss_L=36.0, + K_Pbind=37.2, K_Pdiss=0.93): + ''' + + + Parameters + ---------- + + Returns + ------- + None. + + ''' + pass + @chemicals_user class ModifiedADM1(CompiledProcesses): """ @@ -445,6 +483,25 @@ class ModifiedADM1(CompiledProcesses): KS_IP : float, optional Inorganic phosphorus (nutrient) inhibition coefficient for soluble substrate uptake [M]. The default is 2e-5. + q_aging_H : float, optional + Aging rate constant of X_HFO_H and X_HFO_HP [d^(-1)]. The default is 450.0. + q_aging_L : float, optional + Aging rate constant of X_HFO_L and X_HFO_LP [d^(-1)]. The default is 0.1. + q_Pcoprec : float, optional + Rate constant of P binding and coprecipitation on X_HFO_H [d^(-1)]. + The default is 360. + q_Pbinding : float, optional + Rate constant of P binding on X_HFO_L [d^(-1)]. The default is 0.3. + q_diss_H : float, optional + Dissolution rate constant of X_HFO_HP [d^(-1)]. The default is 36.0. + q_diss_L : float, optional + Dissolution rate constant of X_HFO_HP [d^(-1)]. The default is 36.0. + K_Pbind : float, optional + S_IP half saturation coefficient for binding with X_HFO_H or X_HFO_L + [kg P/m3]. The default is 37.2, i.e., 1.20 kmol P/m3. + K_Pdiss : float, optional + S_IP half inhibition coefficient for dissolution of X_HFO_HP or X_HFO_LP + [kg P/m3]. The default is 0.93, i.e., 0.03 kmol P/m3. KI_h2s_c4 : float, optional H2S half inhibition coefficient for butyrate or valerate uptake [kg COD/m3]. The default is 0.481. @@ -478,7 +535,16 @@ class ModifiedADM1(CompiledProcesses): pH_limits_h2_SRB : 2-tuple, optional Lower and upper limits of pH inhibition for hydrogen uptake by SRB, unitless. The default is (5,6). - + k_cryst : iterable[float], optional + Mineral precipitation rate constants [h^(-1)], following the order of + `ModifiedADM1._precipitates`. The default is + [0.35, 1e-3, 3.0, 1e-3, 2.0, 0.76, 5.0, 1e-3, 1e-3, 1e-3, 1e2, 1e-3, 1e-3]. + n_cryst : iterable[int], optional + The effect orders of mineral precipitation reactions [unitless], following + the order of `ModifiedADM1._precipitates`. The default is + [2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2]. + + Examples -------- ... @@ -518,20 +584,28 @@ class ModifiedADM1(CompiledProcesses): ) _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'K_PP', 'K_so4', 'pH_limits', 'KS_IN', 'KS_IP', 'KI_nh3', 'KIs_h2', 'KIs_h2s' - 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', + 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', + 'k_cryst', 'n_cryst', 'Ksp_base', 'Ksp_dH', 'T_base', 'components', # 'root' ) _acid_base_pairs = ADM1._acid_base_pairs _biogas_IDs = (*ADM1._biogas_IDs, 'S_IS') _biomass_IDs = (*ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') - _precipitate_IDs = ('X_CCM', 'X_ACC', 'X_ACP', 'X_HAP', 'X_DCPD', 'X_OCP', - 'X_struv', 'X_newb', 'X_magn', 'X_kstruv', - 'X_FeS', 'X_Fe3PO42', 'X_AlPO4') + _precipitates = ('X_CCM', 'X_ACC', 'X_ACP', 'X_HAP', 'X_DCPD', 'X_OCP', + 'X_struv', 'X_newb', 'X_magn', 'X_kstruv', + 'X_FeS', 'X_Fe3PO42', 'X_AlPO4') _T_base = 298.15 _K_H_base = [7.8e-4, 1.4e-3, 3.5e-2, 0.105] # biogas species Henry's Law constant [M/bar] _K_H_dH = [-4180, -14240, -19410, -19180] # Heat of reaction of liquid-gas transfer of biogas species [J/mol] + _pKsp_base = [8.48, 8.3, 28.92, 44.333, 18.995, 47.08, + 13.6, 18.175, 7.46, 11.5508, + 2.95, 37.76, 18.2] + _Ksp_dH = [8000, -12000, 54000, 0, 31000, 0, + -22600, -22600, -20000, -22600, + -11000, 5060, 0] + def __new__(cls, components=None, path=None, f_ch_xb=0.275, f_pr_xb=0.275, f_li_xb=0.35, f_xI_xb=0.1, f_fa_li=0.95, f_bu_su=0.1328, f_pro_su=0.2691, f_ac_su=0.4076, @@ -550,13 +624,18 @@ def __new__(cls, components=None, path=None, K_hSRB=5.96e-6, K_aSRB=0.176, K_pSRB=0.088, K_c4SRB=0.1739, K_so4_hSRB=1.04e-4*S_mw, K_so4_aSRB=2e-4*S_mw, K_so4_pSRB=2e-4*S_mw, K_so4_c4SRB=2e-4*S_mw, k_Fe3t2_h2=1e9/Fe_mw, k_Fe3t2_is=1e9/Fe_mw, + q_aging_H=450.0, q_aging_L=0.1, q_Pcoprec=360, q_Pbinding=0.3, q_diss_H=36.0, q_diss_L=36.0, + K_Pbind=37.2, K_Pdiss=0.93, # 1.20 and 0.03 in MATLAB, assuming in kmol-P/m3 ? KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, KS_IP=2e-5, KI_h2s_c4=0.481, KI_h2s_pro=0.481, KI_h2s_ac=0.460, KI_h2s_h2=0.400, KI_h2s_c4SRB=0.520, KI_h2s_pSRB=0.520, KI_h2s_aSRB=0.499, KI_h2s_hSRB=0.499, pH_limits_aa=(4,5.5), pH_limits_ac=(6,7), pH_limits_h2=(5,6), pH_limits_aa_SRB=(6,7), pH_limits_ac_SRB=(6,7), pH_limits_h2_SRB=(5,6), kLa=200, pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], - Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0],**kwargs): + Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0], + k_cryst=[0.35, 1e-3, 3.0, 1e-3, 2.0, 0.76, 5.0, 1e-3, 1e-3, 1e-3, 1e2, 1e-3, 1e-3], + n_cryst=[2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2], + **kwargs): cmps = _load_components(components) @@ -572,7 +651,7 @@ def __new__(cls, components=None, path=None, p.ref_component = 'S_IP' precipitation = [] - for i in cls._precipitate_IDs[:-3]: + for i in cls._precipitates[:-3]: new_p = Process('precipitation_%s' % i.lstrip('X_'), reaction='[?]S_IC + [?]S_IN + [?]S_IP + [?]S_K + [?]S_Mg + [?]S_Ca -> %s' % i, ref_component=i, @@ -644,11 +723,14 @@ def __new__(cls, components=None, path=None, b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2, # original ADM1 q_pha, q_pha, q_pha, q_pha, b_pao, b_pp, b_pha, # P extension k_hSRB, b_hSRB, k_aSRB, b_aSRB, k_pSRB, b_pSRB, k_c4SRB, k_c4SRB, b_c4SRB, # S extension - k_Fe3t2_h2, k_Fe3t2_h2, k_Fe3t2_is, k_Fe3t2_is)) # Fe extension + HFO + k_Fe3t2_h2, k_Fe3t2_h2, k_Fe3t2_is, k_Fe3t2_is, # Fe extension + q_aging_H, q_aging_L, q_Pcoprec, q_Pbinding, # HFO module + q_aging_H, q_aging_L, q_diss_H, q_diss_L)) Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2, # original ADM1 K_A, # P extension - K_hSRB, K_aSRB, K_pSRB, K_c4SRB)) # S extension + K_hSRB, K_aSRB, K_pSRB, K_c4SRB, # S extension + K_Pbind, K_Pdiss)) # HFO module K_so4 = np.array((K_so4_hSRB, K_so4_aSRB, K_so4_pSRB, K_so4_c4SRB)) KIs_h2 = np.array((KI_h2_fa, KI_h2_c4, KI_h2_c4, KI_h2_pro)) @@ -658,18 +740,23 @@ def __new__(cls, components=None, path=None, K_H_dH = np.array(cls._K_H_dH) Ka_base = np.array([10**(-pKa) for pKa in pKa_base]) Ka_dH = np.array(Ka_dH) + k_cryst = np.array(k_cryst) * 24 # converted to d^(-1) + n_cryst = np.array(n_cryst) + Ksp_base = np.array([10**(-pK) for pK in cls.pKsp_base]) + Ksp_dH = np.array(cls.Ksp_dH) # root = TempState() dct = self.__dict__ dct.update(kwargs) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) - # self.set_rate_function(rhos_madm1) - # self.rate_function._params = dict(zip(cls._kinetic_params, - # [ks, Ks, K_PP, K_so4, - # pH_limits, KS_IN*N_mw, KS_IP*P_mw, - # KI_nh3, KIs_h2, KIs_h2s, - # Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, - # cls.T_base, self._components, - # # root, - # ])) + self.set_rate_function(rhos_madm1) + self.rate_function._params = dict(zip(cls._kinetic_params, + [ks, Ks, K_PP, K_so4, + pH_limits, KS_IN*N_mw, KS_IP*P_mw, + KI_nh3, KIs_h2, KIs_h2s, + Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, + k_cryst, n_cryst, Ksp_base, Ksp_dH, + cls.T_base, self._components, + # root, + ])) return self \ No newline at end of file From eaeb620cc4e03cb5766832f6a4b9059b25efc5d5 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 12 Oct 2023 09:52:41 -0400 Subject: [PATCH 13/16] allow `StreamImpactItem` to be created from datasheets --- qsdsan/_impact_item.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/qsdsan/_impact_item.py b/qsdsan/_impact_item.py index c7705908..98fc476a 100644 --- a/qsdsan/_impact_item.py +++ b/qsdsan/_impact_item.py @@ -383,8 +383,12 @@ def get_item(cls, ID): @classmethod def _load_from_df(cls, name, df): if name.lower() == 'info': + if 'kind' not in df.columns: df['kind'] = 'ImpactItem' for num in df.index: - new = cls.__new__(cls) + kind = df.iloc[num].kind.lower() + if kind == 'streamimpactitem': + new = StreamImpactItem.__new__(StreamImpactItem) + else: new = cls.__new__(cls) new.__init__(ID=df.iloc[num].ID, functional_unit=df.iloc[num].functional_unit) else: @@ -409,8 +413,9 @@ def load_from_file(cls, path_or_dict, index_col=None): This Excel should have multiple sheets: - - The "info" sheet should have two columns: "ID" (e.g., Cement) \ - and "functional_unit" (e.g., kg) of different impact items. + - The "info" sheet should have three columns: "ID" (e.g., Cement) \ + "functional_unit" (e.g., kg), and "kind" ("ImpactItem" or "StreamImpactItem") + of different impact items. - The remaining sheets should contain characterization factors of \ impact indicators. From 11bda48a3d215749d476c823cfa9b72c65133dce Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 12 Oct 2023 09:39:01 -0700 Subject: [PATCH 14/16] temperature correction --- qsdsan/processes/_madm1.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index ade77197..c841cdf6 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -263,7 +263,7 @@ def pcm(): Cs = np.empty(38+8) sum_stoichios = np.array([2, 2, 5, 9, 3, 8, 3, 3, 2, 3, 2, 2]) -def rhos_madm1(state_arr, params): +def rhos_madm1(state_arr, params, T_op): ks = params['rate_constants'] Ks = params['half_sat_coeffs'] K_PP = params['K_PP'] @@ -283,8 +283,9 @@ def rhos_madm1(state_arr, params): kLa = params['kLa'] k_cryst = params['k_cryst'] n_cryst = params['n_cryst'] + Kspb = params['Ksp_base'] + Ksp_dH = params['Ksp_dH'] T_base = params['T_base'] - # T_temp = params.pop('T_temp', T_base) Cs[:7] = state_arr[13:20] # original ADM1 processes Cs[7:11] = state_arr[19:23] @@ -327,7 +328,24 @@ def rhos_madm1(state_arr, params): # inhibition factors # ****************** - + unit_conversion = mass2mol_conversion(cmps) + if T_op == T_base: + Ka = Kab + KH = KHb / unit_conversion[[7,8,9,30]] + Ksp = Kspb + else: + T_temp = params.pop('T_op', None) + if T_op == T_temp: + params['T_op'] = T_op + Ka = params['Ka'] + KH = params['KH'] + Ksp = params['Ksp'] + else: + params['T_op'] = T_op + Ka = params['Ka'] = Kab * T_correction_factor(T_base, T_op, Ka_dH) + KH = params['KH'] = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[[7,8,9,30]] + Ksp = params['Ksp'] = Kspb * T_correction_factor(T_base, T_op, Ksp_dH) + S_IN, S_IP = state_arr[[10,11]] I_nutrients = substr_inhibit(S_IN, KS_IN) * substr_inhibit(S_IP, KS_IP) rhos[3:11] *= I_nutrients @@ -360,7 +378,7 @@ def rhos_madm1(state_arr, params): # !!! place holder for PCM (speciation) # ============================================================================= acts = pcm(state_arr) - SIs = np.maximum(1.0, saturation_index(acts)) # should be an array + SIs = np.maximum(1.0, saturation_index(acts, Ksp)) # should be an array rhos[46:59] = k_cryst * state_arr[47:60] * (SIs**(1/sum_stoichios) - 1)**n_cryst #%% modified ADM1 class From 6a549ea950a004e420e6f7e2a8cf47914cc7faa8 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 12 Oct 2023 09:43:38 -0700 Subject: [PATCH 15/16] speed up temperature correction in ADM1 --- qsdsan/processes/_adm1.py | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/qsdsan/processes/_adm1.py b/qsdsan/processes/_adm1.py index cb876244..161f9cfb 100644 --- a/qsdsan/processes/_adm1.py +++ b/qsdsan/processes/_adm1.py @@ -190,6 +190,8 @@ def mass2mol_conversion(cmps): # return np.exp(theta * (T2-T1)) def T_correction_factor(T1, T2, delta_H): + """compute temperature correction factor for equilibrium constants based on + the Van't Holf equation.""" if T1 == T2: return 1 return np.exp(delta_H/(R*100) * (1/T1 - 1/T2)) # R converted to SI @@ -268,10 +270,24 @@ def rhos_adm1(state_arr, params): weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]] T_op = state_arr[-1] + if T_op == T_base: + Ka = Kab + KH = KHb / unit_conversion[7:10] + else: + T_temp = params.pop('T_op', None) + if T_op == T_temp: + params['T_op'] = T_op + Ka = params['Ka'] + KH = params['KH'] + else: + params['T_op'] = T_op + Ka = params['Ka'] = Kab * T_correction_factor(T_base, T_op, Ka_dH) + KH = params['KH'] = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] + biogas_S = state_arr[7:10].copy() biogas_p = R * T_op * state_arr[27:30] - Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) - KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] + # Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) + # KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] rhos[:-3] = ks * Cs Monod = substr_inhibit(substrates, Ks) @@ -280,12 +296,12 @@ def rhos_adm1(state_arr, params): if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu) h = brenth(acid_base_rxn, 1e-14, 1.0, - args=(weak_acids, Kas), + args=(weak_acids, Ka), xtol=1e-12, maxiter=100) # h = 10**(-7.46) - nh3 = Kas[1] * weak_acids[2] / (Kas[1] + h) - co2 = weak_acids[3] - Kas[2] * weak_acids[3] / (Kas[2] + h) + nh3 = Ka[1] * weak_acids[2] / (Ka[1] + h) + co2 = weak_acids[3] - Ka[2] * weak_acids[3] / (Ka[2] + h) biogas_S[-1] = co2 / unit_conversion[9] Iph = Hill_inhibit(h, pH_ULs, pH_LLs) From e71534cf34fd918d64ccdca7b7b1d5f7ab2dacca Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 12 Oct 2023 14:04:06 -0700 Subject: [PATCH 16/16] gas transfer processes --- qsdsan/processes/_madm1.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py index c841cdf6..02ff1b77 100644 --- a/qsdsan/processes/_madm1.py +++ b/qsdsan/processes/_madm1.py @@ -258,6 +258,8 @@ def calc_biogas(): def pcm(): pass +def saturation_index(): + pass rhos = np.zeros(38+8+13+4) # 38 biological + 8 chemical P removal by HFO + 13 MMP + 4 gas transfer Cs = np.empty(38+8) @@ -352,9 +354,9 @@ def rhos_madm1(state_arr, params, T_op): rhos[[25,27,29,31,32]] *= I_nutrients # ============================================================================= -# !!! place holder for pH +# !!! place holder for PCM (speciation) # ============================================================================= - pH, nh3 = calc_pH(state_arr, params) + pH, nh3, co2, acts = pcm(state_arr, params) Is_pH = Hill_inhibit(10**(-pH), pH_ULs, pH_LLs) rhos[3:9] *= Is_pH[0] rhos[9:11] *= Is_pH[1:3] @@ -366,20 +368,25 @@ def rhos_madm1(state_arr, params, T_op): Inh3 = non_compet_inhibit(nh3, KI_nh3) rhos[9] *= Inh3 -# ============================================================================= -# !!! place holder for gas-liquid transfer -# ============================================================================= Z_h2s = calc_biogas() # should be a function of pH, like co2 and nh3 Is_h2s = non_compet_inhibit(Z_h2s, KIs_h2s) rhos[6:11] *= Is_h2s[:5] rhos[[25,27,29,31,32]] *= Is_h2s[5:] -# ============================================================================= -# !!! place holder for PCM (speciation) -# ============================================================================= - acts = pcm(state_arr) + # multiple mineral precipitation + # ****************************** SIs = np.maximum(1.0, saturation_index(acts, Ksp)) # should be an array rhos[46:59] = k_cryst * state_arr[47:60] * (SIs**(1/sum_stoichios) - 1)**n_cryst + + # gas transfer + # ************ + biogas_S = state_arr[[7,8,9,30]].copy() + biogas_S[2] = co2 / unit_conversion[9] + biogas_S[3] = Z_h2s / unit_conversion[30] + biogas_p = R * T_op * state_arr[63:67] + rhos[-4:] = kLa * (biogas_S - KH * biogas_p) + + return rhos #%% modified ADM1 class _load_components = settings.get_default_chemicals