enable general pH control in AnaerobicCSTR

QSD-Group · Jul 29, 2024 · fbfb59a · fbfb59a
1 parent 112b7df
commit fbfb59a
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 4 deletions.
diff --git a/qsdsan/processes/_adm1.py b/qsdsan/processes/_adm1.py
@@ -629,6 +629,7 @@ def __new__(cls, components=None, path=None, N_xc=2.686e-3, N_I=4.286e-3, N_aa=7
                                                K_H_base, K_H_dH, kLa,
                                                T_base, self._components, root]))
 
+        dct['flex_rhos'] = _rhos_adm1
         dct['solve_pH'] = solve_pH
         dct['dydt_Sh2_AD'] = dydt_Sh2_AD
         dct['grad_dydt_Sh2_AD'] = grad_dydt_Sh2_AD

diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py
@@ -463,6 +463,7 @@ def __new__(cls, components=None, path=None,
                                                T_base, self._components, root, 
                                                #!!! new parameter
                                                KS_IP*P_mw]))
+        dct['flex_rhos'] = _rhos_adm1_p_extension
         dct['solve_pH'] = solve_pH
         dct['dydt_Sh2_AD'] = dydt_Sh2_AD
         dct['grad_dydt_Sh2_AD'] = grad_dydt_Sh2_AD
@@ -1041,6 +1042,7 @@ def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in):
             Q = state_arr[45]
             return -Q/V_liq + np.dot(grad_rhos, stoichio[[5,6,7,8,10]]) + kLa*stoichio[-3]
 
+        dct['flex_rhos'] = _rhos_adm1p
         dct['solve_pH'] = adm1p_solve_pH
         dct['dydt_Sh2_AD'] = dydt_Sh2_AD
         dct['grad_dydt_Sh2_AD'] = grad_dydt_Sh2_AD

diff --git a/qsdsan/sanunits/_anaerobic_reactor.py b/qsdsan/sanunits/_anaerobic_reactor.py
@@ -270,7 +270,8 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None,
                  T=308.15, headspace_P=1.013, external_P=1.013, 
                  pipe_resistance=5.0e4, fixed_headspace_P=False,
                  retain_cmps=(), fraction_retain=0.95,
-                 isdynamic=True, exogenous_vars=(), **kwargs):
+                 isdynamic=True, exogenous_vars=(), 
+                 pH_ctrl=None, **kwargs):
         if len(exogenous_vars) == 0:
             exogenous_vars = (EDV('T', function=lambda t: T), )
         super().__init__(ID=ID, ins=ins, outs=outs, thermo=thermo,
@@ -293,6 +294,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None,
         self.fixed_headspace_P = fixed_headspace_P
         self._f_retain = np.array([fraction_retain if cmp.ID in retain_cmps \
                                    else 0 for cmp in self.components])
+        self.pH_ctrl = pH_ctrl
         self._mixed = WasteStream()
         self._tempstate = {}
 
@@ -523,10 +525,10 @@ 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.algebraic_h2)
+            self._compile_ODE(self.algebraic_h2, self.pH_ctrl)
         return self._ODE
 
-    def _compile_ODE(self, algebraic_h2=True):
+    def _compile_ODE(self, algebraic_h2=True, pH_ctrl=None):
         if self._model is None:
             CSTR._compile_ODE(self)
         else:
@@ -535,7 +537,11 @@ def _compile_ODE(self, algebraic_h2=True):
             _state = self._state
             _dstate = self._dstate
             _update_dstate = self._update_dstate
-            _f_rhos = self.model.rate_function
+            if pH_ctrl:
+                _params = self.model.rate_function.params
+                _f_rhos = lambda state_arr: self.model.flex_rhos(state_arr, _params, h=10**(-pH_ctrl))
+            else:
+                _f_rhos = self.model.rate_function
             _f_param = self.model.params_eval
             _M_stoichio = self.model.stoichio_eval
             n_cmps = len(cmps)