enable DEA solving with AnaerobicCSTR

qsdsan/processes/_adm1.py

@@ -26,7 +26,9 @@
            'non_compet_inhibit', 'substr_inhibit',
            'T_correction_factor', 
            'pH_inhibit', 'Hill_inhibit', 
-           'rhos_adm1')
+           'rhos_adm1', 
+           'solve_pH',
+           'dydt_Sh2_AD', 'grad_dydt_Sh2_AD')
 
 _path = ospath.join(data_path, 'process_data/_adm1.tsv')
 _load_components = settings.get_default_chemicals
@@ -325,7 +327,7 @@ def _rhos_adm1(state_arr, params, h=None):
     rhos[-3:] = kLa * (biogas_S - KH * biogas_p)
     return rhos
 
-def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, S_h2_in, V_liq):
+def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in):
     state_arr[7] = S_h2
     Q = state_arr[30]
     rxn = _rhos_adm1(state_arr, params, h=h)
@@ -334,10 +336,11 @@ def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, S_h2_in, V_liq):
 
 grad_rhos = np.zeros(5)
 X_bio = np.zeros(5)
-def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq):
+def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in):
     state_arr[7] = S_h2
     ks = params['rate_constants'][[6,7,8,9,11]]
-    Ks = params['half_sat_coeffs'][6:10]
+    Ks = params['half_sat_coeffs'][2:6]
+    K_h2 = params['half_sat_coeffs'][7]
     pH_ULs = params['pH_ULs']
     pH_LLs = params['pH_LLs']
     KS_IN = params['KS_IN']
@@ -349,14 +352,13 @@ def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq):
     S_va, S_bu, S_IN = state_arr[[3,4,10]]
     Iph = Hill_inhibit(h, pH_ULs, pH_LLs)[[2,3,4,5,7]]
     Iin = substr_inhibit(S_IN, KS_IN)
-    gIh2 = grad_non_compet_inhibit(S_h2, KIs_h2)
+    grad_Ih2 = grad_non_compet_inhibit(S_h2, KIs_h2)
 
     grad_rhos[:] = ks * X_bio * Iph * Iin
-    grad_rhos[:4] *= substr_inhibit(substrates, Ks) * gIh2
+    grad_rhos[:-1] *= substr_inhibit(substrates, Ks) * grad_Ih2
     if S_va > 0: rhos[1] *= 1/(1+S_bu/S_va)
     if S_bu > 0: rhos[2] *= 1/(1+S_va/S_bu)
 
-    K_h2 = params['half_sat_coeffs'][11]
     grad_rhos[-1] *= grad_substr_inhibit(S_h2, K_h2)
     stoichio = f_stoichio(state_arr)
 

qsdsan/sanunits/_anaerobic_reactor.py

@@ -18,10 +18,13 @@
 from math import ceil, pi
 from biosteam import Stream
 from .. import SanUnit, Construction, WasteStream
-from ..processes import Decay
+from ..processes import Decay, T_correction_factor, solve_pH, \
+    dydt_Sh2_AD, grad_dydt_Sh2_AD
 from ..sanunits import HXutility, WWTpump, CSTR
 from ..utils import ospath, load_data, data_path, auom, \
     calculate_excavation_volume, ExogenousDynamicVariable as EDV
+from scipy.optimize import newton
+
 __all__ = (
     'AnaerobicBaffledReactor',
     'AnaerobicCSTR',
@@ -261,6 +264,7 @@ class AnaerobicCSTR(CSTR):
     _ins_size_is_fixed = False
     _outs_size_is_fixed = False
     _R = 8.3145e-2 # Universal gas constant, [bar/M/K]
+    algebraic_h2 = True
 
     def __init__(self, ID='', ins=None, outs=(), thermo=None,
                  init_with='WasteStream', V_liq=3400, V_gas=300, model=None,  
@@ -519,15 +523,16 @@ def f_q_gas_var_P_headspace(self, rhoTs, S_gas, T):
     @property
     def ODE(self):
         if self._ODE is None:
-            self._compile_ODE()
+            self._compile_ODE(self.algebraic_h2)
         return self._ODE
 
-    def _compile_ODE(self):
+    def _compile_ODE(self, algebraic_h2=True):
         if self._model is None:
             CSTR._compile_ODE(self)
         else:
             cmps = self.components
             f_rtn = self._f_retain
+            _state = self._state
             _dstate = self._dstate
             _update_dstate = self._update_dstate
             _f_rhos = self.model.rate_function
@@ -544,37 +549,64 @@ def _compile_ODE(self):
                 f_qgas = self.f_q_gas_fixed_P_headspace
             else:
                 f_qgas = self.f_q_gas_var_P_headspace
-
             if self.model._dyn_params:
                 def M_stoichio(state_arr):
                     _f_param(state_arr)
                     return self.model.stoichio_eval().T
             else:
                 _M_stoichio = self.model.stoichio_eval().T
                 M_stoichio = lambda state_arr: _M_stoichio
+
+            h2_idx = cmps.index('S_h2')
+            if algebraic_h2:
+                params = self.model.rate_function.params
+                if self.model._dyn_params:
+                    def h2_stoichio(state_arr):
+                        return M_stoichio(state_arr)[h2_idx]
+                else:
+                    _h2_stoichio = _M_stoichio[h2_idx]
+                    h2_stoichio = lambda state_arr: _h2_stoichio
+                unit_conversion = cmps.i_mass / cmps.chem_MW
+                def solve_h2(QC, S_in, T):
+                    Ka = params['Ka_base'] * T_correction_factor(params['T_base'], T, params['Ka_dH'])
+                    h, nh3, co2 = solve_pH(QC, Ka, unit_conversion)
+                    S_h2_0 = QC[h2_idx]
+                    S_h2_in = S_in[h2_idx]
+                    S_h2 = newton(dydt_Sh2_AD, S_h2_0, grad_dydt_Sh2_AD,
+                                  args=(QC, h, params, h2_stoichio, V_liq, S_h2_in), 
+                                  )
+                    return S_h2
+                def update_h2_dstate(dstate):
+                    dstate[h2_idx] = 0.
+            else:
+                solve_h2 = lambda QC, S_ins, T: QC[h2_idx]
+                def update_h2_dstate(dstate):
+                    pass
             def dy_dt(t, QC_ins, QC, dQC_ins):
                 QC[QC < 2.2e-16] = 0.
-                S_liq = QC[:n_cmps]
-                S_gas = QC[n_cmps: (n_cmps+n_gas)]
-                #!!! Volume change due to temperature difference accounted for 
-                # in _run and _init_state
                 Q_ins = QC_ins[:, -1]
                 S_ins = QC_ins[:, :-1] * 1e-3  # mg/L to kg/m3
                 Q = sum(Q_ins)
+                S_in = Q_ins @ S_ins / Q
                 if hasexo:
                     exo_vars = f_exovars(t)
                     QC = np.append(QC, exo_vars)
                     T = exo_vars[0]
                 else: T = self.T
+                QC[h2_idx] = _state[h2_idx] = solve_h2(QC, S_in, T)
                 rhos =_f_rhos(QC)
+                S_liq = QC[:n_cmps]
+                S_gas = QC[n_cmps: (n_cmps+n_gas)]
                 _dstate[:n_cmps] = (Q_ins @ S_ins - Q*S_liq*(1-f_rtn))/V_liq \
                     + np.dot(M_stoichio(QC), rhos)
                 q_gas = f_qgas(rhos[-3:], S_gas, T)
                 _dstate[n_cmps: (n_cmps+n_gas)] = - q_gas*S_gas/V_gas \
                     + rhos[-3:] * V_liq/V_gas * gas_mass2mol_conversion
                 # _dstate[-1] = dQC_ins[0,-1]
                 _dstate[-1] = 0.
+                update_h2_dstate(_dstate)
                 _update_dstate()
+
             self._ODE = dy_dt
 
     def get_retained_mass(self, biomass_IDs):

qsdsan/sanunits/_sludge_treatment.py

@@ -352,8 +352,6 @@ def AE(self):
 
     def _compile_AE(self):
 
-        # This function is run multiple times during dynamic simulation 
-
         _state = self._state
         _dstate = self._dstate
         _update_state = self._update_state