modified ADM1 complete component set

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

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__,

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 <[email protected]>
+
+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
+
+#%%

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