From 1e1f966bec1fee221bb046d09175565583e6412f Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 12:04:17 -0700 Subject: [PATCH 01/14] merge changes from main --- docs/source/conf.py | 2 +- qsdsan/processes/__init__.py | 6 +- qsdsan/processes/_adm1_p_extension.py | 46 +- qsdsan/processes/_aeration.py | 21 + qsdsan/processes/_asm2d.py | 46 +- .../sanunits/_suspended_growth_bioreactor.py | 511 +++++++++--------- tests/test_junctions.py | 497 +++++++++++++++++ 7 files changed, 868 insertions(+), 261 deletions(-) create mode 100644 tests/test_junctions.py diff --git a/docs/source/conf.py b/docs/source/conf.py index 7ea627c7..cc735083 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -31,7 +31,7 @@ # built documents. # # The short X.Y version. -version = '1.3.1' +version = '1.4.0' # The full version, including alpha/beta/rc tags. release = version diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index 5aacbc5a..ba9beecb 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -63,7 +63,7 @@ def __init__(self): from ._asm2d import * from ._adm1 import * from ._adm1_p_extension import * -from ._madm1 import * +# from ._madm1 import * from ._decay import * from ._kinetic_reaction import * from ._pm2 import * @@ -74,7 +74,7 @@ def __init__(self): _asm2d, _adm1, _adm1_p_extension, - _madm1, + # _madm1, _decay, _kinetic_reaction, _pm2 @@ -86,7 +86,7 @@ def __init__(self): *_asm2d.__all__, *_adm1.__all__, *_adm1_p_extension.__all__, - *_madm1.__all__, + # *_madm1.__all__, *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py index 1c1a81cb..aa4d1818 100644 --- a/qsdsan/processes/_adm1_p_extension.py +++ b/qsdsan/processes/_adm1_p_extension.py @@ -94,15 +94,6 @@ def acid_base_rxn(h_ion, weak_acids_tot, Kas): nh3, hpo4, hco3, ac, pro, bu, va = Kas[1:] * weak_acids_tot[4:] / (Kas[1:] + h_ion) return S_cat + S_K + 2*S_Mg + h_ion + (S_IN - nh3) - S_an - oh_ion - hco3 - ac - pro - bu - va - 2*hpo4 - (S_IP - hpo4) -# The function 'fprime_abr' is not used in the code -def fprime_abr(h_ion, weak_acids_tot, Kas): - # S_cat, S_K, S_Mg, S_an, S_IN, S_IP = weak_acids_tot[:6] - Kw = Kas[0] - doh_ion = - Kw / h_ion ** 2 - dnh3, dhpo4, dhco3, dac, dpro, dbu, dva = - Kas[1:] * weak_acids_tot[4:] / (Kas[1:] + h_ion)**2 - return 1 + (-dnh3) - doh_ion - dhco3 - dac - dpro - dbu - dva - dhpo4 - - rhos = np.zeros(28) # 28 kinetic processes (25 as defined in modified ADM1 + 3 for gases) Cs = np.empty(25) # 25 processes as defined in modified ADM1 @@ -239,7 +230,6 @@ def _rhos_adm1_p_extension(state_arr, params, h=None): rhos[9] *= Inh3 rhos[-3:] = kLa * (biogas_S - KH * biogas_p) - # print(rhos) return rhos def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): @@ -355,6 +345,42 @@ class ADM1_p_extension(ADM1): >>> adm1_p.show() ADM1_p_extension([hydrolysis_carbs, hydrolysis_proteins, hydrolysis_lipids, uptake_sugars, uptake_amino_acids, uptake_LCFA, uptake_valerate, uptake_butyrate, uptake_propionate, uptake_acetate, uptake_h2, decay_Xsu, decay_Xaa, decay_Xfa, decay_Xc4, decay_Xpro, decay_Xac, decay_Xh2, storage_Sva_in_XPHA, storage_Sbu_in_XPHA, storage_Spro_in_XPHA, storage_Sac_in_XPHA, lysis_XPAO, lysis_XPP, lysis_XPHA, h2_transfer, ch4_transfer, IC_transfer]) + >>> import numpy as np + >>> state_arr = np.ones(cmps.size + len(adm1_p._biogas_IDs) + 2) # liquid-phase concentrations, gas-phase concentrations, liquid flowrate, and temperature + >>> state_arr[-1] = 273.15+35 # Temperature + >>> rhos = adm1_p.rate_function(state_arr) # reaction rate for each process + >>> for i,j in zip(adm1_p.IDs, rhos): + ... print(f'{i}{(40-len(i))*" "}{j:.3g}') + hydrolysis_carbs 10 + hydrolysis_proteins 10 + hydrolysis_lipids 10 + uptake_sugars 20 + uptake_amino_acids 38.4 + uptake_LCFA 2.14e-05 + uptake_valerate 8.32e-05 + uptake_butyrate 8.32e-05 + uptake_propionate 4.13e-05 + uptake_acetate 1.93 + uptake_h2 34.9 + decay_Xsu 0.02 + decay_Xaa 0.02 + decay_Xfa 0.02 + decay_Xc4 0.02 + decay_Xpro 0.02 + decay_Xac 0.02 + decay_Xh2 0.02 + storage_Sva_in_XPHA 0.747 + storage_Sbu_in_XPHA 0.747 + storage_Spro_in_XPHA 0.747 + storage_Sac_in_XPHA 0.747 + lysis_XPAO 0.2 + lysis_XPP 0.2 + lysis_XPHA 0.2 + h2_transfer 139 + ch4_transfer -181 + IC_transfer -1.66e+03 + + References ---------- [1] Batstone, D. J.; Keller, J.; Angelidaki, I.; Kalyuzhnyi, S. V; diff --git a/qsdsan/processes/_aeration.py b/qsdsan/processes/_aeration.py index d5a0ffec..9da6d678 100644 --- a/qsdsan/processes/_aeration.py +++ b/qsdsan/processes/_aeration.py @@ -80,6 +80,25 @@ class DiffusedAeration(Process): [parameters] KLa: 240 DOsat: 8 [dynamic parameters] + + >>> aer2 = pc.DiffusedAeration('aer2', 'S_O', KLa_20=100, V=1000, d_submergence=3.7) + >>> aer2.show() + Process: aer2 + [stoichiometry] S_O: 1 + [reference] S_O + [rate equation] KLa*(DOsat - S_O) + [parameters] KLa: 60 + DOsat: 9.87 + [dynamic parameters] + + >>> aer2.Q_air # doctest: +ELLIPSIS + 12470.65... + >>> round(aer2.SOTR / 1000) + 1039 + + >>> aer3 = pc.DiffusedAeration('aer3', 'S_O', Q_air=7600, V=1000, d_submergence=3.7) + >>> aer3.kLa # doctest: +ELLIPSIS + 36.56... """ @@ -357,5 +376,7 @@ def DOsat(self): return self._DOsat @DOsat.setter def DOsat(self, DOsat): + if DOsat is not None: + self._DOsat_s20 = DOsat / (self.tau * self._beta * self.Omega * self.delta) self._DOsat = DOsat or self._calc_DOsat() self.set_parameters(DOsat = self._DOsat) \ No newline at end of file diff --git a/qsdsan/processes/_asm2d.py b/qsdsan/processes/_asm2d.py index d12e90c5..c2a99b7e 100644 --- a/qsdsan/processes/_asm2d.py +++ b/qsdsan/processes/_asm2d.py @@ -73,8 +73,6 @@ def create_asm2d_cmps(set_thermo=True): return cmps_asm2d -# create_asm2d_cmps() - def create_masm2d_cmps(set_thermo=True): c2d = create_asm2d_cmps(False) ion_kwargs = dict(particle_size='Soluble', @@ -433,7 +431,7 @@ class ASM2d(CompiledProcesses): Examples -------- - >>> from qsdsan import processes as pc, set_thermo + >>> from qsdsan import processes as pc >>> cmps = pc.create_asm2d_cmps() >>> asm2d = pc.ASM2d() >>> asm2d.show() @@ -714,6 +712,8 @@ class mASM2d(CompiledProcesses): electron_acceptor_dependent_decay : bool, optional Whether biomass decay kinetics is dependent on concentrations of electron acceptors. The default is True. + pH_ctrl : float or None, optional + Whether to fix pH at a specific value or solve for pH (`None`). The default is 7.0. k_h : float, optional Hydrolysis rate constant, in [d^(-1)]. The default is 3.0. eta_NO3_Hl : float, optional @@ -766,6 +766,46 @@ class mASM2d(CompiledProcesses): >>> asm.show() mASM2d([aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, storage_PHA, aero_storage_PP, anox_storage_PP, PAO_aero_growth_PHA, PAO_anox_growth, PAO_lysis, PP_lysis, PHA_lysis, auto_aero_growth, auto_lysis, CaCO3_precipitation_dissolution, struvite_precipitation_dissolution, newberyite_precipitation_dissolution, ACP_precipitation_dissolution, MgCO3_precipitation_dissolution, AlPO4_precipitation_dissolution, FePO4_precipitation_dissolution]) + >>> # Calculate process rate given state variable values and fixed pH. + >>> import numpy as np + >>> state_arr = np.ones(len(cmps)) + >>> rhos = asm.rate_function(state_arr) # reaction rate for each process + >>> for i,j in zip(asm.IDs, rhos): + ... print(f'{i}{(40-len(i))*" "}{j:.3g}') + aero_hydrolysis 2.27 + anox_hydrolysis 0.182 + anae_hydrolysis 0.0606 + hetero_growth_S_F 0.471 + hetero_growth_S_A 0.471 + denitri_S_F 0.0503 + denitri_S_A 0.0503 + ferment 0.0333 + hetero_lysis 0.356 + storage_PHA 0.594 + aero_storage_PP 1.06 + anox_storage_PP 0.0851 + PAO_aero_growth_PHA 0.778 + PAO_anox_growth 0.0622 + PAO_lysis 0.174 + PP_lysis 0.174 + PHA_lysis 0.174 + auto_aero_growth 0.33 + auto_lysis 0.111 + CaCO3_precipitation_dissolution 0 + struvite_precipitation_dissolution 0 + newberyite_precipitation_dissolution 0 + ACP_precipitation_dissolution 0 + MgCO3_precipitation_dissolution 0 + AlPO4_precipitation_dissolution 1.82e-11 + FePO4_precipitation_dissolution 1.82e-11 + + >>> # Estimate pH given state variable values. + >>> Ka = asm.rate_function.params['Ka'] + >>> unit_conversion = asm.rate_function.params['mass2mol'] + >>> h_ion = asm.solve_pH(state_arr, Ka, unit_conversion) + >>> pH = -np.log10(h_ion) + >>> print(f'{pH:.2f}') + 8.40 References ---------- diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 5c305a59..4b8cc3e8 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -10,20 +10,19 @@ for license details. ''' -from .. import SanUnit, WasteStream, Process, Processes, CompiledProcesses -from ._clarifier import _settling_flux +from .. import SanUnit, WasteStream, Process, CompiledProcesses +# from ._clarifier import _settling_flux from ..sanunits import dydt_cstr -from sympy import symbols, lambdify, Matrix -from scipy.integrate import solve_ivp +# from scipy.integrate import solve_ivp from warnings import warn -from math import floor, ceil +# from math import floor, ceil import numpy as np import pandas as pd -from numba import njit +# from numba import njit __all__ = ('CSTR', 'BatchExperiment', - 'SBR', + # 'SBR', 'PFR', ) @@ -437,7 +436,40 @@ def _cost(self): #%% class BatchExperiment(SanUnit): + ''' + A batch reactor in experimental settings. + + Parameters + ---------- + model : :class:`CompiledProcesses`, optional + Process model that describes the dynamics of state variables. + The `state` of the batch reactor is entirely determined by the + stoichiometry and rate function in this model. + + Examples + -------- + >>> import qsdsan.sanunits as su, qsdsan.processes as pc + >>> cmps = pc.create_asm1_cmps() + >>> asm1 = pc.ASM1() + >>> BE = su.BatchExperiment('BE', model=asm1) + >>> BE.set_init_conc(S_S=20, X_BH=500, S_O=8, S_ND=3, S_ALK=84) + >>> BE.simulate(t_span=(0,10), method='BDF') + >>> for k,v in BE.state.items(): + ... if v != 0: + ... print(f'{k}{" "*(7-len(k))}{v:.2f}') + S_S 0.93 + X_S 446.16 + X_BH 25.66 + X_P 39.25 + S_O -0.00 + S_NO -0.00 + S_NH 2.12 + S_ND 0.00 + X_ND 36.47 + S_ALK 85.82 + S_N2 0.00 + ''' _N_ins = 0 _N_outs = 0 # _ins_size_is_fixed = True @@ -445,16 +477,7 @@ class BatchExperiment(SanUnit): def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', model=None, isdynamic=True, exogenous_vars=(), **kwargs): - ''' - A batch reactor in experimental settings. - - Parameters - ---------- - model : :class:`CompiledProcesses`, optional - Process model that describes the dynamics of state variables. - The `state` of the batch reactor is entirely determined by the - stoichiometry and rate function in this model. - ''' + SanUnit.__init__(self, ID, None, (), thermo, init_with, isdynamic=isdynamic, exogenous_vars=exogenous_vars) self._model = model @@ -530,233 +553,233 @@ def dy_dt(t, QC_ins, QC, dQC_ins): #TODO: add functions for convenient model calibration #%% NOT READY -class SBR(SanUnit): - ''' - Sequential batch reactors operated in parallel. The number of reactors is - determined by operation cycle and influent flowrate. [1]_ - - Parameters - ---------- - ID : str - ID for the reactors. The default is ''. - ins : :class:`WasteStream` - Influent to the reactor. Expected number of influent is 1. - outs : :class:`WasteStream` - Treated effluent and wasted sludge. - surface_area : float, optional - Surface area of the reactor bottom, in [m^2]. The reactor is assumed - to be cylinder. The default is 1500. - height : float, optional - Height of the reactor, in [m]. The default is 4. - operation_cycle : iterable of float, optional - Operation cycle of the SBR, time for each stage specified in [h]. There - are 7 stages: 1 - fill, 2 - fill, 3 - mix, 4 - mix, 5 - settle, 6 - decant, - 7 - desludge. The first 4 stages are modeled as a biological reactor. - The 5th stage is modeled as a 1D N-layer settler. The last 2 stages are - assumed inactive. The default is (0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1). - aeration : iterable of float and/or :class:`Process`, optional - Aeration settings for the first 4 stages of the operation cycle. Either - specify a targeted dissolved oxygen concentration in [mg O2/L] or provide - a :class:`Process` object to represent aeration, or None for no aeration. - The default is (None, None, None, 2.0). - DO_ID : str, optional - The :class:`Component` ID for dissolved oxygen, only relevant when the - reactor is aerated. The default is 'S_O2'. - suspended_growth_model : :class:`Processes`, optional - The suspended growth biokinetic model. The default is None. - N_layer : int, optional - The number of layers to model settling. The default is 10. - pumped_flow : float, optional - Designed effluent flowrate, in [m^3/d]. The default is None. - underflow : float, optional - Designed wasted activated sludge flowrate, in [m^3/d]. The default is None. - X_threshold : float, optional - Threshold suspended solid concentration, in [g/m^3]. The default is 3000. - v_max : float, optional - Maximum theoretical (i.e. Vesilind) settling velocity, in [m/d]. The - default is 474. - v_max_practical : float, optional - Maximum practical settling velocity, in [m/d]. The default is 250. - rh : float, optional - Hindered zone settling parameter in the double-exponential settling velocity - function, in [m^3/g]. The default is 5.76e-4. - rp : float, optional - Flocculant zone settling parameter in the double-exponential settling velocity - function, in [m^3/g]. The default is 2.86e-3. - fns : float, optional - Non-settleable fraction of the suspended solids, dimensionless. Must be within - [0, 1]. The default is 2.28e-3. - cache_state : bool, optional - Whether to store volume and composition of retained sludge in the tank from - most recent run. The default is True. - - References - ---------- - .. [1] Takács, I.; Patry, G. G.; Nolasco, D. A Dynamic Model of the Clarification - -Thickening Process. Water Res. 1991, 25 (10), 1263–1271. - https://doi.org/10.1016/0043-1354(91)90066-Y. - - ''' - - _N_ins = 1 - _N_outs = 2 - - def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', - surface_area=1500, height=4, - operation_cycle=(0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1), - aeration=(None, None, None, 2.0), DO_ID='S_O2', - suspended_growth_model=None, N_layer=10, - pumped_flow=None, underflow=None, - X_threshold=3000, v_max=474, v_max_practical=250, - rh=5.76e-4, rp=2.86e-3, fns=2.28e-3, - cache_state=True, **kwargs): - SanUnit.__init__(self, ID, ins, outs, thermo, init_with) - - self._V = surface_area * height - self._A = surface_area - self._h = height - self._operation_cycle = operation_cycle - self._aeration = aeration - self._DO_ID = DO_ID - self._model = suspended_growth_model - self._N_layer = N_layer - self._Q_e = pumped_flow - self._Q_WAS = underflow - self._X_t = X_threshold - self._v_max = v_max - self._v_max_p = v_max_practical - self._rh = rh - self._rp = rp - self._fns = fns - self._cache_state = cache_state - - for attr, value in kwargs.items(): - setattr(self, attr, value) - self._init_Vas = None - self._init_Cas = None - self._dynamic_composition = None - - - @property - def operation_cycle(self): - return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2', 'settle', 'decant', 'desludge'), - self._operation_cycle)) - @property - def total_cycle_time(self): - return sum(self._operation_cycle) - - @property - def aeration(self): - return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2'), - self._aeration[:4])) - - @property - def C_t(self): - if self._dynamic_composition: - return pd.DataFrame(self._dynamic_composition, - columns = ['Time[d]'] + list(self.components.IDs)) - else: return None - - def _run(self, cache_state=True): - if self._model is None: - raise RuntimeError(f'{self.ID} was initialized without a suspended growth model.') - else: - isa = isinstance - inf = self.ins[0] - Q_in = inf.get_total_flow('m3/d') - eff, sludge = self.outs - eff.copy_like(inf) - sludge.copy_like(inf) - C_in = inf.mass / inf.F_vol * 1e3 # concentrations in g/m3 - cmps = self.components - C = list(symbols(cmps.IDs)) - if self._init_Vas is not None: - V_0 = self._init_Vas - C_0 = self._init_Cas - else: - V_0 = 0 - C_0 = C_in - n = self._N_layer - if self._aeration.count(None) == len(self._aeration): - Vmax = self._V - hj = self._h/n - else: - Vmax = self._V*0.75 - hj = self._h*0.75/n - - # ********fill and mix/aerate stages*********** - T_fill = (Vmax - V_0)/Q_in # maximum total fill time in day - T = [t/24 for t in self._operation_cycle] # operation cycle in day - if T_fill <= T[0]: - schedule = [T_fill, T[0]-T_fill] + T[1:4] - aer = [self._aeration[0], self._aeration[0]] + list(self._aeration[1:4]) - fill = [True] + [False]*4 - V_total = Vmax - elif T_fill <= T[0]+T[1]: - schedule = [T[0], T_fill-T[0], T[0]+T[1]-T_fill] + T[2:4] - aer = list(self._aeration[:2]) + [self._aeration[1]] + list(self._aeration[2:4]) - fill = [True]*2 + [False]*3 - V_total = Vmax - else: - schedule = T[:4] - aer = list(self._aeration[:4]) - fill = [True]*2 + [False]*2 - V_total = Q_in*(T[0]+T[1])+V_0 - hj = V_total/self._V*self._h/n - - for i in range(1, len(schedule)): - if fill[-i] == fill[-i-1] and aer[-i] == aer[-i-1]: - schedule[-i-1] += schedule[-i] - schedule[-i] = 0 - - t_arr = np.array([]) - y_mat = np.ndarray([]) - for i in range(len(schedule)): - if schedule[i] > 0: - dC_dt, J_func = self._compile_dC_dt(V_0, Q_in, C_in, C, fill[i], aer[i]) - if isa(aer[i], (float, int)): C_0[cmps.index(self._DO_ID)] = aer[i] - sol = solve_ivp(dC_dt, (0, schedule[i]), C_0, method='BDF', jac=J_func) - C_0 = sol.y.transpose()[-1] - V_0 += Q_in * schedule[i] * fill[i] - t_arr = np.concatenate((t_arr, sol.t + t_arr[-1])) - y_mat = np.hstack((y_mat, sol.y)) - self._dynamic_composition = np.vstack((t_arr, y_mat)).transpose() - - # *********settle, decant, desludge********** - eff.set_flow(C_0*eff.F_vol, 'g/hr', self.components.IDs) - X_0 = eff.get_TSS() - X_min = X_0 * self._fns - T_settle = T[4] - def dX_dt(t, X): - VX = [_settling_flux(x, self._v_max, self._v_max_p, X_min, self._rh, self._rp) for x in X] - J = [VX[j] if X[j+1] <= self._X_t else min(VX[j], VX[j+1]) for j in range(n-1)] - settle_out = np.array(J + [0]) - settle_in = np.array([0] + J) - dXdt = (settle_in - settle_out)/hj - return dXdt - sol = solve_ivp(dX_dt, (0, T_settle), np.ones(n)*X_0) - X = sol.y.transpose()[-1] - - V_eff = min(T[5]*self._Q_e, V_total*(n-1)/n) - n_eff = V_eff/V_total - w_eff = np.array([1]*floor(n_eff)+[n_eff-floor(n_eff)]) - X_eff = np.average(X[:ceil(n_eff)], weights=w_eff) - eff_mass_flow = (X_eff/X_0*cmps.x + (1-cmps.x))*C_0*V_eff/T[5] - eff.set_flow(eff_mass_flow, 'g/d', cmps.IDs) - - V_was = min(T[6]*self._Q_WAS, V_total-V_eff) - X_as = (V_total*X_0 - V_eff*X_eff) / (V_total-V_eff) - C_as = (X_as/X_0*cmps.x + (1-cmps.x))*C_0 - was_mass_flow = C_as*V_was/T[6] - sludge.set_flow(was_mass_flow, 'g/d', cmps.IDs) - - if self._cache_state: - self._init_Vas = V_total - V_eff - V_was - self._init_Cas = C_as - - - def _design(self): - pass +# class SBR(SanUnit): +# ''' +# Sequential batch reactors operated in parallel. The number of reactors is +# determined by operation cycle and influent flowrate. [1]_ + +# Parameters +# ---------- +# ID : str +# ID for the reactors. The default is ''. +# ins : :class:`WasteStream` +# Influent to the reactor. Expected number of influent is 1. +# outs : :class:`WasteStream` +# Treated effluent and wasted sludge. +# surface_area : float, optional +# Surface area of the reactor bottom, in [m^2]. The reactor is assumed +# to be cylinder. The default is 1500. +# height : float, optional +# Height of the reactor, in [m]. The default is 4. +# operation_cycle : iterable of float, optional +# Operation cycle of the SBR, time for each stage specified in [h]. There +# are 7 stages: 1 - fill, 2 - fill, 3 - mix, 4 - mix, 5 - settle, 6 - decant, +# 7 - desludge. The first 4 stages are modeled as a biological reactor. +# The 5th stage is modeled as a 1D N-layer settler. The last 2 stages are +# assumed inactive. The default is (0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1). +# aeration : iterable of float and/or :class:`Process`, optional +# Aeration settings for the first 4 stages of the operation cycle. Either +# specify a targeted dissolved oxygen concentration in [mg O2/L] or provide +# a :class:`Process` object to represent aeration, or None for no aeration. +# The default is (None, None, None, 2.0). +# DO_ID : str, optional +# The :class:`Component` ID for dissolved oxygen, only relevant when the +# reactor is aerated. The default is 'S_O2'. +# suspended_growth_model : :class:`Processes`, optional +# The suspended growth biokinetic model. The default is None. +# N_layer : int, optional +# The number of layers to model settling. The default is 10. +# pumped_flow : float, optional +# Designed effluent flowrate, in [m^3/d]. The default is None. +# underflow : float, optional +# Designed wasted activated sludge flowrate, in [m^3/d]. The default is None. +# X_threshold : float, optional +# Threshold suspended solid concentration, in [g/m^3]. The default is 3000. +# v_max : float, optional +# Maximum theoretical (i.e. Vesilind) settling velocity, in [m/d]. The +# default is 474. +# v_max_practical : float, optional +# Maximum practical settling velocity, in [m/d]. The default is 250. +# rh : float, optional +# Hindered zone settling parameter in the double-exponential settling velocity +# function, in [m^3/g]. The default is 5.76e-4. +# rp : float, optional +# Flocculant zone settling parameter in the double-exponential settling velocity +# function, in [m^3/g]. The default is 2.86e-3. +# fns : float, optional +# Non-settleable fraction of the suspended solids, dimensionless. Must be within +# [0, 1]. The default is 2.28e-3. +# cache_state : bool, optional +# Whether to store volume and composition of retained sludge in the tank from +# most recent run. The default is True. + +# References +# ---------- +# .. [1] Takács, I.; Patry, G. G.; Nolasco, D. A Dynamic Model of the Clarification +# -Thickening Process. Water Res. 1991, 25 (10), 1263–1271. +# https://doi.org/10.1016/0043-1354(91)90066-Y. + +# ''' + +# _N_ins = 1 +# _N_outs = 2 + +# def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', +# surface_area=1500, height=4, +# operation_cycle=(0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1), +# aeration=(None, None, None, 2.0), DO_ID='S_O2', +# suspended_growth_model=None, N_layer=10, +# pumped_flow=None, underflow=None, +# X_threshold=3000, v_max=474, v_max_practical=250, +# rh=5.76e-4, rp=2.86e-3, fns=2.28e-3, +# cache_state=True, **kwargs): +# SanUnit.__init__(self, ID, ins, outs, thermo, init_with) + +# self._V = surface_area * height +# self._A = surface_area +# self._h = height +# self._operation_cycle = operation_cycle +# self._aeration = aeration +# self._DO_ID = DO_ID +# self._model = suspended_growth_model +# self._N_layer = N_layer +# self._Q_e = pumped_flow +# self._Q_WAS = underflow +# self._X_t = X_threshold +# self._v_max = v_max +# self._v_max_p = v_max_practical +# self._rh = rh +# self._rp = rp +# self._fns = fns +# self._cache_state = cache_state + +# for attr, value in kwargs.items(): +# setattr(self, attr, value) +# self._init_Vas = None +# self._init_Cas = None +# self._dynamic_composition = None + + +# @property +# def operation_cycle(self): +# return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2', 'settle', 'decant', 'desludge'), +# self._operation_cycle)) +# @property +# def total_cycle_time(self): +# return sum(self._operation_cycle) + +# @property +# def aeration(self): +# return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2'), +# self._aeration[:4])) + +# @property +# def C_t(self): +# if self._dynamic_composition: +# return pd.DataFrame(self._dynamic_composition, +# columns = ['Time[d]'] + list(self.components.IDs)) +# else: return None + +# def _run(self, cache_state=True): +# if self._model is None: +# raise RuntimeError(f'{self.ID} was initialized without a suspended growth model.') +# else: +# isa = isinstance +# inf = self.ins[0] +# Q_in = inf.get_total_flow('m3/d') +# eff, sludge = self.outs +# eff.copy_like(inf) +# sludge.copy_like(inf) +# C_in = inf.mass / inf.F_vol * 1e3 # concentrations in g/m3 +# cmps = self.components +# C = list(symbols(cmps.IDs)) +# if self._init_Vas is not None: +# V_0 = self._init_Vas +# C_0 = self._init_Cas +# else: +# V_0 = 0 +# C_0 = C_in +# n = self._N_layer +# if self._aeration.count(None) == len(self._aeration): +# Vmax = self._V +# hj = self._h/n +# else: +# Vmax = self._V*0.75 +# hj = self._h*0.75/n + +# # ********fill and mix/aerate stages*********** +# T_fill = (Vmax - V_0)/Q_in # maximum total fill time in day +# T = [t/24 for t in self._operation_cycle] # operation cycle in day +# if T_fill <= T[0]: +# schedule = [T_fill, T[0]-T_fill] + T[1:4] +# aer = [self._aeration[0], self._aeration[0]] + list(self._aeration[1:4]) +# fill = [True] + [False]*4 +# V_total = Vmax +# elif T_fill <= T[0]+T[1]: +# schedule = [T[0], T_fill-T[0], T[0]+T[1]-T_fill] + T[2:4] +# aer = list(self._aeration[:2]) + [self._aeration[1]] + list(self._aeration[2:4]) +# fill = [True]*2 + [False]*3 +# V_total = Vmax +# else: +# schedule = T[:4] +# aer = list(self._aeration[:4]) +# fill = [True]*2 + [False]*2 +# V_total = Q_in*(T[0]+T[1])+V_0 +# hj = V_total/self._V*self._h/n + +# for i in range(1, len(schedule)): +# if fill[-i] == fill[-i-1] and aer[-i] == aer[-i-1]: +# schedule[-i-1] += schedule[-i] +# schedule[-i] = 0 + +# t_arr = np.array([]) +# y_mat = np.ndarray([]) +# for i in range(len(schedule)): +# if schedule[i] > 0: +# dC_dt, J_func = self._compile_dC_dt(V_0, Q_in, C_in, C, fill[i], aer[i]) +# if isa(aer[i], (float, int)): C_0[cmps.index(self._DO_ID)] = aer[i] +# sol = solve_ivp(dC_dt, (0, schedule[i]), C_0, method='BDF', jac=J_func) +# C_0 = sol.y.transpose()[-1] +# V_0 += Q_in * schedule[i] * fill[i] +# t_arr = np.concatenate((t_arr, sol.t + t_arr[-1])) +# y_mat = np.hstack((y_mat, sol.y)) +# self._dynamic_composition = np.vstack((t_arr, y_mat)).transpose() + +# # *********settle, decant, desludge********** +# eff.set_flow(C_0*eff.F_vol, 'g/hr', self.components.IDs) +# X_0 = eff.get_TSS() +# X_min = X_0 * self._fns +# T_settle = T[4] +# def dX_dt(t, X): +# VX = [_settling_flux(x, self._v_max, self._v_max_p, X_min, self._rh, self._rp) for x in X] +# J = [VX[j] if X[j+1] <= self._X_t else min(VX[j], VX[j+1]) for j in range(n-1)] +# settle_out = np.array(J + [0]) +# settle_in = np.array([0] + J) +# dXdt = (settle_in - settle_out)/hj +# return dXdt +# sol = solve_ivp(dX_dt, (0, T_settle), np.ones(n)*X_0) +# X = sol.y.transpose()[-1] + +# V_eff = min(T[5]*self._Q_e, V_total*(n-1)/n) +# n_eff = V_eff/V_total +# w_eff = np.array([1]*floor(n_eff)+[n_eff-floor(n_eff)]) +# X_eff = np.average(X[:ceil(n_eff)], weights=w_eff) +# eff_mass_flow = (X_eff/X_0*cmps.x + (1-cmps.x))*C_0*V_eff/T[5] +# eff.set_flow(eff_mass_flow, 'g/d', cmps.IDs) + +# V_was = min(T[6]*self._Q_WAS, V_total-V_eff) +# X_as = (V_total*X_0 - V_eff*X_eff) / (V_total-V_eff) +# C_as = (X_as/X_0*cmps.x + (1-cmps.x))*C_0 +# was_mass_flow = C_as*V_was/T[6] +# sludge.set_flow(was_mass_flow, 'g/d', cmps.IDs) + +# if self._cache_state: +# self._init_Vas = V_total - V_eff - V_was +# self._init_Cas = C_as + + +# def _design(self): +# pass # def _compile_dC_dt(self, V0, Qin, Cin, C, fill, aer): # isa = isinstance diff --git a/tests/test_junctions.py b/tests/test_junctions.py new file mode 100644 index 00000000..3edad88b --- /dev/null +++ b/tests/test_junctions.py @@ -0,0 +1,497 @@ +# -*- coding: utf-8 -*- +''' +EXPOsan: Exposition of sanitation and resource recovery systems + +This module is developed by: + + Joy Zhang + + Yalin Li + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt +for license details. + +Reference: +.. [1] Alex, J.; Benedetti, L.; Copp, J. B.; Gernaey, K. V.; Jeppsson, U.; + Nopens, I.; Pons, M. N.; Rosen, C.; Steyer, J. P.; Vanrolleghem, P. A. + Benchmark Simulation Model No. 2 (BSM2). + http://iwa-mia.org/wp-content/uploads/2022/09/TR3_BSM_TG_Tech_Report_no_3_BSM2_General_Description.pdf. +.. [2] 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 + +''' +#%% + +def test_adm1_junctions(): + + import qsdsan as qs, numpy as np + from numpy.testing import assert_allclose as ac + from qsdsan import ( + processes as pc, + sanunits as su, + WasteStream, + ) + from qsdsan.utils import ospath, load_data + from exposan.bsm2 import data_path + + matlab_preAD_adm = { + 'S_su': 0.0, # monosacharides (kg COD/m3) + 'S_aa': 0.04388, # amino acids (kg COD/m3) + 'S_fa': 0.0, # long chain fatty acids (LCFA) (kg COD/m3) + 'S_va': 0.0, # total valerate (kg COD/m3) + 'S_bu': 0.0, # total butyrate (kg COD/m3) + 'S_pro': 0.0, # total propionate (kg COD/m3) + 'S_ac': 0.0, # total acetate (kg COD/m3) + 'S_h2': 0.0, # hydrogen gas (kg COD/m3) + 'S_ch4': 0.0, # methane gas (kg COD/m3) + 'S_IC': 0.0079326*12, # inorganic carbon (kmole C/m3 -> kg C/m3) 0.0951912 + 'S_IN': 0.0019721*14, # inorganic nitrogen (kmole N/m3 -> kg N/m3) 0.0276094 + 'S_I': 0.028067, # soluble inerts (kg COD/m3) + 'X_c': 0.0, # composites (kg COD/m3) + 'X_ch': 3.7236, # carbohydrates (kg COD/m3) + 'X_pr': 15.9235, # proteins (kg COD/m3) + 'X_li': 8.047, # lipids (kg COD/m3) + 'X_su': 0.0, # sugar degraders (kg COD/m3) + 'X_aa': 0.0, # amino acid degraders (kg COD/m3) + 'X_fa': 0.0, # LCFA degraders (kg COD/m3) + 'X_c4': 0.0, # valerate and butyrate degraders (kg COD/m3) + 'X_pro': 0.0, # propionate degraders (kg COD/m3) + 'X_ac': 0.0, # acetate degraders (kg COD/m3) + 'X_h2': 0.0, # hydrogen degraders (kg COD/m3) + 'X_I': 17.0106, # particulate inerts (kg COD/m3) + 'S_cat': 0.0, # cations (base) (kmole/m3) + 'S_an': 0.0052101, # anions (acid) (kmole/m3) + # 'Q': 178.4674, # Flow rate (m3/d) + } + + matlab_postAD_adm = { + 'S_su': 0.012394, + 'S_aa': 0.0055432, + 'S_fa': 0.10741, + 'S_va': 0.012333, + 'S_bu': 0.014003, + 'S_pro': 0.017584, + 'S_ac': 0.089315, + 'S_h2': 2.5055e-07, + 'S_ch4': 0.05549, + 'S_IC': 0.095149*12, + 'S_IN': 1.3226, + 'S_I': 0.13087, + 'X_c': 0.10792, + 'X_ch': 0.020517, + 'X_pr': 0.08422, + 'X_li': 0.043629, + 'X_su': 0.31222, + 'X_aa': 0.93167, + 'X_fa': 0.33839, + 'X_c4': 0.33577, + 'X_pro': 0.10112, + 'X_ac': 0.67724, + 'X_h2': 0.28484, + 'X_I': 17.2162, + 'S_cat': 0., #-4.0789e-34, + 'S_an': 0.0052101 + } + + matlab_postAD_asm = { + 'S_I': 130.867, # soluble inert organic matter, mg COD/l + 'S_S': 258.5789, # readily biodegradable substrate, mg COD/l + 'X_I': 17216.2434, # particulate inert organic matter, mg COD/l + 'X_S': 2611.4843, # slowly biodegradable substrate, mg COD/l + 'X_BH': 0.0, # active heterotrophic biomass, mg COD/l + 'X_BA': 0.0, # active autotrophic biomass, mg COD/l + 'X_P': 626.0652, # particulate products arising from biomass decay, mg COD/l + 'S_O': 0.0, # dissolved O2, mg -COD/l + 'S_NO': 0.0, # nitrate and nitrite nitrogen, mg N/L + 'S_NH': 1442.7882, # ammonium, mg N/L + 'S_ND': 0.54323, # soluble biodegradable organic nitrogen + 'X_ND': 100.8668, # particulate biodegradable organic nitrogen, mg N/l + 'S_ALK': 97.8459*12, # alkalinity, assumed to be HCO3-, 97.8459, mol HCO3/m3 -> g C/m3 + 'S_N2': 0.0, # dissolved O2 + # 'Q': 178.4674, # Flow rate, m3/d + } + + + adm1init = load_data(ospath.join(data_path, 'adm1init.csv'), index_col=0).to_dict('index') + asm1_default_parameters = dict( + mu_H = 4.0, + K_S = 10.0, + K_OH = 0.2, + K_NO = 0.5, + b_H = 0.3, + mu_A = 0.5, + K_NH = 1.0, + K_OA = 0.4, + b_A = 0.05, + eta_g = 0.8, + k_a = 0.05, + k_h = 3.0, + K_X = 0.1, + eta_h = 0.8, + Y_H = 0.67, + Y_A = 0.24, + f_P = 0.08, + i_XB = 0.08, + i_XP = 0.06, + fr_SS_COD = 0.75 + ) + + T = 273.15 + 35 + cmps_asm1 = pc.create_asm1_cmps() + asm1 = pc.ASM1(components=cmps_asm1, **asm1_default_parameters) + preAD_asm = WasteStream('preAD_asm', T=T) + preAD_asm.set_flow_by_concentration( + flow_tot=178.4674, + concentrations=dict( + S_I = 28.0665, + S_S = 48.9526, + X_I = 10361.7101, + X_S = 20375.0176, + X_BH = 10210.0698, + X_BA = 553.2808, + X_P = 3204.6601, + S_O = 0.25225, + S_NO = 1.6871, + S_NH = 28.9098, + S_ND = 4.6834, + X_ND = 906.0933, + S_ALK = 7.1549*12 + ), + units=('m3/d', 'mg/L') + ) + thermo_asm1 = qs.get_thermo() + cmps_adm1 = pc.create_adm1_cmps() + adm1 = pc.ADM1() + cmps_adm1.X_I.i_N = cmps_asm1.X_I.i_N # slight difference + cmps_adm1.refresh_constants() + thermo_adm1 = qs.get_thermo() + + J1 = su.ASMtoADM('J1', upstream=preAD_asm, downstream='preAD_adm', + thermo=thermo_adm1, isdynamic=True, adm1_model=adm1,#) + T=T, pH=7.2631) + AD1 = su.AnaerobicCSTR('AD1', ins=J1-0, outs=('biogas', 'postAD_adm'), + isdynamic=True, V_liq=3400, V_gas=300, T=T, + model=adm1,) + AD1.set_init_conc(**adm1init['AD1']) + # Switch back to ASM1 components + J2 = su.ADMtoASM('J2', upstream=AD1-1, downstream='postAD_asm', + thermo=thermo_asm1, isdynamic=True, adm1_model=adm1) + J2.bio_to_xs = 0.79 + qs.set_thermo(thermo_asm1) + + sys = qs.System(path=(J1, AD1, J2)) + sys.simulate(state_reset_hook='reset_cache', t_span=(0, 200), method='BDF') + fs = sys.flowsheet.stream + + for ws in sys.streams: + ws.state[ws.state < 2.2e-16] = 0 + + ac(cmps_adm1.kwarray(matlab_preAD_adm)[:-1]*1e3, fs.preAD_adm.state[:-2], rtol=1e-4) + ac(cmps_adm1.kwarray(matlab_postAD_adm)[:-1]*1e3, fs.postAD_adm.state[:-2], rtol=1e-2) + ac(cmps_asm1.kwarray(matlab_postAD_asm)[:-1], fs.postAD_asm.state[:-2], rtol=1e-3) + + h2 = cmps_adm1.S_h2 + ch4 = cmps_adm1.S_ch4 + co2 = cmps_adm1.S_IC + assert np.isclose(AD1.state['S_h2_gas'] * h2.chem_MW / h2.i_mass, 1.1032e-5, rtol=1e-3) + assert np.isclose(AD1.state['S_ch4_gas'] * ch4.chem_MW / ch4.i_mass, 1.6535, rtol=1e-2) + assert np.isclose(AD1.state['S_IC_gas'], 0.01354, rtol=1e-2) + assert np.isclose(AD1.outs[1].pH, 7.2631, rtol=1e-3) + + assert np.isclose(fs.biogas.imass['S_h2']*24 * h2.i_mass, 0.0035541, rtol=1e-2) + assert np.isclose(fs.biogas.imass['S_ch4']*24 * ch4.i_mass, 1065.3523, rtol=1e-2) + assert np.isclose(fs.biogas.imass['S_IC']*24 * co2.i_mass, 1535.4118, rtol=1e-2) + + sys.flowsheet.clear() + +#%% +def test_adm1p_junctions(): + import numpy as np + from numpy.testing import assert_allclose as ac + from chemicals.elements import molecular_weight as get_mw + from qsdsan import sanunits as su, processes as pc, WasteStream, System, get_thermo + # from qsdsan.utils import load_data, ospath, time_printer + # from exposan.bsm2 import data_path + + Q = 190 # influent flowrate [m3/d] + HRT = 20 + V_liq = Q*HRT + V_gas = 0.088*V_liq + Temp = 273.15+35 # temperature [K] + C_mw = get_mw({'C':1}) + N_mw = get_mw({'N':1}) + P_mw = get_mw({'P':1}) + struv_mw = get_mw(dict(Mg=1, N=1, H=4, P=1, O=4)) + # adm1init = load_data(ospath.join(data_path, 'adm1init.csv'), index_col=0).to_dict('index') + + # Table 1.1 [mg/L], Flores-Alsina et al., 2016. Appendix + inf_asm2d = dict( + S_O2=0, + S_F=26.44, + S_A=17.66, + S_I=27.23, + S_NH4=18.58, + S_N2=5.07, + S_NO3=0.02, + S_PO4=4.69, + S_IC=78.99, + X_I=10964.41, + X_S=19084.76, + X_H=9479.39, + X_PAO=3862.20, + X_PP=450.87, + X_PHA=24.64, + X_AUT=333.79, + S_K=19.79, + S_Mg=189.87, + S_Na=70, + S_Cl=1035, + S_Ca=300, + ) + + # Table 1.3 [kg/m3] + inf_adm1p = dict( + S_su=0.018, + S_aa=0.008, + S_ac=0.018, + S_IC=0.021*C_mw, + S_IN=0.036*N_mw, + S_IP=0.006*P_mw, + S_I=0.027, + X_ch=8.020, + X_pr=8.481, + X_li=11.416, + X_I=11.946, + X_PHA=0.025, + X_PP=0.015*P_mw, + X_PAO=3.862, + S_K=0.001*39, + S_Mg=0.008*24.3, + S_Ca=0.007*40, + S_Na=0.003*23, + S_Cl=0.029*35.5, + # S_N2=0.0004*14 + ) + + # [kmol/m3] + _inf_adm1p = dict( + S_IC=0.021, + S_IN=0.036, + S_IP=0.006, + X_PP=0.015, + S_K=0.001, + S_Mg=0.008, + S_Ca=0.007, + S_Na=0.003, + S_Cl=0.029, + ) + + # Table 1.4 [kg/m3] + out_adm1p = dict( + S_su=0.013, + S_aa=0.006, + S_fa=0.116, + S_va=0.012, + S_bu=0.016, + S_pro=0.019, + S_ac=0.055, + S_h2=2.65e-7, + S_ch4=0.052, + S_IC=0.059*C_mw, + S_IN=0.080*N_mw, + S_IP=0.007*P_mw, + S_I=0.027, + X_ch=1.441, + X_pr=1.513, + X_li=2.025, + X_I=12.345, + X_PHA=0.252, + X_PP=8.05e-6*P_mw, + # X_biomass=3.600, + X_su=3.600, + S_K=0.005*39, + S_Mg=0.001*24.3, + S_Ca=0.001*40, + X_ACP=0.002*310.176722, + X_struv=0.011*245.406502, + S_Na=0.003*23, + S_Cl=0.029*35.5, + # S_N2=0.0004*14 + ) + + # _out_adm1p = dict( + # S_IC=0.059, + # S_IN=0.080, + # S_IP=0.007, + # X_PP=8.05e-6, + # S_K=0.005, + # S_Mg=0.001, + # S_Ca=0.001, + # X_ACP=0.002, + # X_struv=0.011, + # S_Na=0.003, + # S_Cl=0.029, + # ) + + # Table 1.5 [mg/L] + out_asm2d = dict( + S_NH4=1291.68, + S_PO4=298.09, + S_F=134.43, + S_A=353.82, + S_I=27.23, + S_IC=885.27, + S_K=208.84, + S_Mg=28.29, + X_I=12704.93, + X_S=8218.94, + S_Na=70, + S_Cl=1035, + S_Ca=20.45, + X_ACP=722.17, + X_struv=1578.52*245.406502/struv_mw + ) + + # [mmol/L] + _out_asm2d = dict( + S_NH4=1291.68/N_mw, + S_PO4=298.09/P_mw, + S_IC=885.27/C_mw, + S_K=208.84/39, + S_Mg=28.29/24.3, + S_Na=70/23, + S_Cl=1035/35.5, + S_Ca=20.45/40, + X_ACP=722.17/310.176722, + X_struv=1578.52/struv_mw + ) + + default_init_conds = { + 'S_su': 0.014*1e3, + 'S_aa': 0.0062*1e3, + 'S_fa': 0.126*1e3, + 'S_va': 0.0129*1e3, + 'S_bu': 0.0168*1e3, + 'S_pro': 0.0204*1e3, + 'S_ac': 0.0588*1e3, + 'S_h2': 2.8309e-7*1e3, + 'S_ch4': 0.0544*1e3, + 'S_IC': 0.089*12*1e3, + 'S_IN': 0.0663*14*1e3, + 'S_IP': 0.028*31*1e3, + 'S_I': 0.1309*1e3, + 'X_ch': 1.302*1e3, + 'X_pr': 1.3613*1e3, + 'X_li': 1.8127*1e3, + 'X_su': 0.5146*1e3, + 'X_aa': 0.4017*1e3, + 'X_fa': 0.3749*1e3, + 'X_c4': 0.1596*1e3, + 'X_pro': 0.0896*1e3, + 'X_ac': 0.5006*1e3, + 'X_h2': 0.258*1e3, + 'X_I': 12.9232*1e3, + 'X_PHA': 0.6697*1e3, + 'X_PAO': 0.9154*1e3, + 'S_K': 0.0129*1e3, + 'S_Mg': 0.0001*1e3, + 'S_Ca': 2e-4*1e3, + 'X_struv':0.0161*1e3, + 'X_ACP': 9e-4*1e3, + 'X_FePO4': 0.001*1e3, + 'S_Na': 0.061*1e3, + 'S_Cl': 0.0126*1e3 + } + + cmps_asm = pc.create_masm2d_cmps() + inf_asm = WasteStream('inf_asm', T=Temp) + inf_asm.set_flow_by_concentration( + flow_tot=Q, + concentrations=inf_asm2d, + units=('m3/d', 'mg/L') + ) + alt_eff_asm = WasteStream('alt_eff_asm', T=Temp) + alt_eff_asm.set_flow_by_concentration( + flow_tot=Q, + concentrations=out_asm2d, + units=('m3/d', 'mg/L') + ) + asm = pc.mASM2d() + thermo_asm = get_thermo() + cmps_adm = pc.create_adm1p_cmps() + alt_inf_adm = WasteStream('alt_inf_adm', T=Temp) + alt_inf_adm.set_flow_by_concentration( + flow_tot=Q, + concentrations=inf_adm1p, + units=('m3/d', 'kg/m3') + ) + alt_eff_adm = WasteStream('alt_eff_adm', T=Temp) + alt_eff_adm.set_flow_by_concentration( + flow_tot=Q, + concentrations=out_adm1p, + units=('m3/d', 'kg/m3') + ) + adm = pc.ADM1p( + f_bu_su=0.1328, f_pro_su=0.2691, f_ac_su=0.4076, + q_ch_hyd=0.3, q_pr_hyd=0.3, q_li_hyd=0.3, + ) + thermo_adm = get_thermo() + + J1 = su.mASM2dtoADM1p( + 'J1', upstream=inf_asm, downstream='inf_adm', + thermo=thermo_adm, isdynamic=True, + adm1_model=adm, asm2d_model=asm + ) + J1.xs_to_li = 0.6 + AD = su.AnaerobicCSTR( + 'AD', + ins=alt_inf_adm, + # ins=J1-0, + outs=('biogas', 'eff_adm'), isdynamic=True, + V_liq=V_liq, V_gas=V_gas, T=Temp, model=adm + ) + AD.algebraic_h2 = False + AD.set_init_conc(**default_init_conds) + J2 = su.ADM1ptomASM2d( + 'J2', + upstream=alt_eff_adm, + # upstream=AD-1, + downstream='eff_asm', thermo=thermo_asm, isdynamic=True, + adm1_model=adm, asm2d_model=asm + ) + + sys = System(path=(J1, AD, J2)) + sys.simulate(state_reset_hook='reset_cache', t_span=(0, 200), method='BDF') + s = sys.flowsheet.stream + + ########## mASM2d to ADM1p ########### + mass2mol = cmps_adm.i_mass / cmps_adm.chem_MW + idx = cmps_adm.indices(_inf_adm1p.keys()) + _molar = np.round(s.inf_adm.conc[idx] * mass2mol[idx] * 1e-3, 3) + ac(_molar, np.array([v for v in _inf_adm1p.values()])) + ac(np.delete(s.inf_adm.conc, idx)[:-1], # exclude water + np.delete(s.alt_inf_adm.conc, idx)[:-1], + atol=1.0) + + ########## !!! ADM1p skip for now ########## + + ######### ADM1p to mASM2d ########### + mass2mol = cmps_asm.i_mass / cmps_asm.chem_MW + idx = cmps_asm.indices(_out_asm2d.keys()) + _molar = s.eff_asm.conc[idx] * mass2mol[idx] + ac(_molar, np.array([v for v in _out_asm2d.values()]), atol=1.0) + ac(np.delete(s.eff_asm.conc, idx)[:-1], # exclude water + np.delete(s.alt_eff_asm.conc, idx)[:-1], + atol=1.0) + + sys.flowsheet.clear() + +#%% + +if __name__ == '__main__': + test_adm1_junctions() + test_adm1p_junctions() From 59f0837a10b50224749792a42060a86286eee1c4 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 24 Oct 2024 15:22:49 -0700 Subject: [PATCH 02/14] `ExcretionmASM2d` excreta i.f.o. `mASM2d` components --- qsdsan/sanunits/_excretion.py | 146 +++++++++++++++++++++++++++++++++- 1 file changed, 143 insertions(+), 3 deletions(-) diff --git a/qsdsan/sanunits/_excretion.py b/qsdsan/sanunits/_excretion.py index 07f8eaa4..22bf3064 100644 --- a/qsdsan/sanunits/_excretion.py +++ b/qsdsan/sanunits/_excretion.py @@ -5,7 +5,10 @@ QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems This module is developed by: + Yalin Li + + 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 @@ -16,8 +19,11 @@ from .. import SanUnit from ..utils import ospath, load_data, data_path +from warnings import warn +# from scipy.linalg import solve as la_solve +import numpy as np -__all__ = ('Excretion',) +__all__ = ('Excretion', 'ExcretionmASM2d') excretion_path = ospath.join(data_path, 'sanunit_data/_excretion.tsv') @@ -44,7 +50,7 @@ class Excretion(SanUnit): [1] Trimmer et al., Navigating Multidimensional Social–Ecological System Trade-Offs across Sanitation Alternatives in an Urban Informal Settlement. Environ. Sci. Technol. 2020, 54 (19), 12641–12653. - https://doi.org/10.1021/acs.est.0c03296. + https://doi.org/10.1021/acs.est.0c03296 ''' _N_ins = 0 @@ -58,6 +64,8 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream data = load_data(path=excretion_path) for para in data.index: value = float(data.loc[para]['expected']) + # value = float(data.loc[para]['low']) + # value = float(data.loc[para]['high']) setattr(self, '_'+para, value) del data @@ -317,4 +325,136 @@ def waste_ratio(self): return self._waste_ratio @waste_ratio.setter def waste_ratio(self, i): - self._waste_ratio = i \ No newline at end of file + self._waste_ratio = i + + +#%% + +class ExcretionmASM2d(Excretion): + + def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', + waste_ratio=0, **kwargs): + super().__init__(ID, ins, outs, thermo, init_with, waste_ratio, **kwargs) + isdyn = kwargs.pop('isdynamic', False) + if isdyn: self._init_dynamic() + + def _run(self): + ur, fec = self.outs + ur.empty() + fec.empty() + cmps = ur.components + sf_iN = cmps.S_F.i_N + xs_iN = cmps.X_S.i_N + xb_iN = cmps.X_H.i_N + sxi_iN = cmps.S_I.i_N + i_mass = cmps.i_mass + i_P = cmps.i_P + hco3_imass = cmps.S_IC.i_mass + + not_wasted = 1 - self.waste_ratio + factor = 24 * 1e3 # from g/cap/d to kg/hr(/cap) + e_cal = self.e_cal / 24 * not_wasted # kcal/cap/d --> kcal/cap/hr + ur_exc = self.ur_exc / factor + fec_exc = self.fec_exc / factor + + # 14 kJ/g COD, the average lower heating value of excreta + tot_COD = e_cal*self.e_exc*4.184/14/1e3 # in kg COD/hr + fec_COD = tot_COD*self.e_fec + ur_COD = tot_COD - fec_COD + + tot_N = (self.p_veg+self.p_anim)*self.N_prot/factor \ + * self.N_exc*not_wasted + ur_N = tot_N*self.N_ur + fec_N = tot_N - ur_N + + tot_P = (self.p_veg*self.P_prot_v+self.p_anim*self.P_prot_a)/factor \ + * self.P_exc*not_wasted + ur_P = tot_P*self.P_ur + fec_P = tot_P - ur_P + + # breakpoint() + ur.imass['S_NH4'] = ur_nh4 = ur_N * self.N_ur_NH3 + req_sf_cod = (ur_N - ur_nh4) / sf_iN + if req_sf_cod <= ur_COD: + ur.imass['S_F'] = sf = req_sf_cod + ur.imass['S_A'] = ur_COD - sf # contains no N or P + else: + req_si_cod = (ur_N - ur_nh4) / sxi_iN + if req_si_cod <= ur_COD: + ur.imass['S_F'] = sf = (sxi_iN * ur_COD - (ur_N - ur_nh4))/(sxi_iN - sf_iN) + ur.imass['S_I'] = ur_COD - sf + else: + ur.imass['S_F'] = sf = ur_COD + ur_other_n = ur_N - ur_nh4 - sf * sf_iN + warn(f"Excess non-NH3 nitrogen cannot be accounted for by organics " + f"in urine: {ur_other_n} kg/hr. Added to NH3-N.") + ur.imass['S_NH4'] += ur_other_n # debatable, has negative COD # raise warning/error + + ur.imass['S_PO4'] = ur_P - sum(ur.mass * i_P) + ur.imass['S_K'] = e_cal/1e3 * self.K_cal/1e3 * self.K_exc*self.K_ur + ur.imass['S_Mg'] = self.Mg_ur / factor + ur.imass['S_Ca'] = self.Ca_ur / factor + + ur.imass['H2O'] = self.ur_moi * ur_exc + ur_others = ur_exc - sum(ur.mass * i_mass) + ur.imass['S_IC'] = ur_others * 0.34 / hco3_imass + ur.imass['S_Na'] = ur_others * 0.35 + ur.imass['S_Cl'] = ur_others * 0.31 + + fec.imass['S_NH4'] = fec_nh4 = fec_N * self.N_fec_NH3 + req_xs_cod = (fec_N - fec_nh4) / xs_iN + if req_xs_cod <= fec_COD: + fec.imass['X_S'] = xs = req_xs_cod + fec.imass['S_A'] = fec_COD - xs + else: + req_xi_cod = (fec_N - fec_nh4) / sxi_iN + if req_xi_cod <= fec_COD: + fec.imass['X_S'] = xs = (sxi_iN * fec_COD - (fec_N - fec_nh4))/(sxi_iN - xs_iN) + fec.imass['X_I'] = fec_COD - xs + else: + req_xb_cod = (fec_N - fec_nh4) / xb_iN + if req_xb_cod <= fec_COD: + fec.imass['X_S'] = xs = (xb_iN * fec_COD - (fec_N - fec_nh4))/(xb_iN - xs_iN) + fec.imass['X_H'] = fec_COD - xs + else: + fec.imass['X_S'] = xs = fec_COD + fec_other_n = fec_N - fec_nh4 - xs * xs_iN + warn(f"Excess non-NH3 nitrogen cannot be accounted for by organics " + f"in feces: {fec_other_n} kg/hr. Added to NH3-N.") + fec.imass['S_NH4'] += fec_other_n # debatable, has negative COD + + fec.imass['S_PO4'] = fec_P - sum(fec.mass * i_P) + fec.imass['S_K'] = (1-self.K_ur)/self.K_ur * ur.imass['S_K'] + fec.imass['S_Mg'] = self.Mg_fec / factor + fec.imass['S_Ca'] = self.Ca_fec / factor + fec.imass['H2O'] = self.fec_moi * fec_exc + + fec_others = fec_exc - sum(fec.mass * i_mass) + fec.imass['S_IC'] = fec_others * 0.34 / hco3_imass + fec.imass['S_Na'] = fec_others * 0.35 + fec.imass['S_Cl'] = fec_others * 0.31 + + + @property + def AE(self): + if self._AE is None: + self._compile_AE() + return self._AE + + def _compile_AE(self): + def yt(t, QC_ins, dQC_ins): + pass + self._AE = yt + + def _init_state(self): + ur, fec = self.outs + self._state = np.append(ur.mass, fec.mass) + for ws in self.outs: + ws.state = np.append(ws.conc, ws.F_vol * 24) + ws.dstate = np.zeros_like(ws.state) + + def _update_state(self): + pass + + def _update_dstate(self): + pass \ No newline at end of file From 90768b20530938f83d2c4bbc7faecfaff81ae650 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 19 Nov 2024 13:28:31 -0800 Subject: [PATCH 03/14] created `CompletelyMixedMBR` --- qsdsan/sanunits/_membrane_bioreactor.py | 184 +++++++++++++++++++++++- 1 file changed, 181 insertions(+), 3 deletions(-) diff --git a/qsdsan/sanunits/_membrane_bioreactor.py b/qsdsan/sanunits/_membrane_bioreactor.py index 1bd90c4f..65e07212 100644 --- a/qsdsan/sanunits/_membrane_bioreactor.py +++ b/qsdsan/sanunits/_membrane_bioreactor.py @@ -5,8 +5,12 @@ QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems This module is developed by: + Yalin Li + 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. @@ -16,7 +20,10 @@ from biosteam import Stream from biosteam.exceptions import DesignError from . import HXutility, WWTpump, InternalCirculationRx -from .. import SanStream, SanUnit +from .. import SanStream, WasteStream, SanUnit, Process, CompiledProcesses +from ..sanunits import CSTR, PFR, dydt_cstr +from warnings import warn +import numpy as np, pandas as pd from ..equipments import Blower from ..utils import ( auom, @@ -26,8 +33,11 @@ calculate_excavation_volume, ) -__all__ = ('AnMBR',) +__all__ = ('AnMBR', + 'CompletelyMixedMBR', + 'PlugFlowMBR',) +#%% degassing = SanStream.degassing _ft_to_m = auom('ft').conversion_factor('m') @@ -1338,4 +1348,172 @@ def organic_rm(self): Qi, Qe = self._inf.F_vol, self.outs[1].F_vol Si = compute_stream_COD(self._inf, 'kg/m3') Se = compute_stream_COD(self.outs[1], 'kg/m3') - return 1 - Qe*Se/(Qi*Si) \ No newline at end of file + return 1 - Qe*Se/(Qi*Si) + +#%% +from numba import njit +@njit(cache=True) +def dydt_mbr(QC_ins, QC, V, Qp, f_rtn, xarr, _dstate): + Q_ins = QC_ins[:, -1] + C_ins = QC_ins[:, :-1] + Q = sum(Q_ins) + C = QC[:-1] + _dstate[-1] = 0 + _dstate[:-1] = (Q_ins @ C_ins - (Q*(1-xarr) + (Qp+(Q-Qp)*(1-f_rtn))*xarr)*C)/V + +class CompletelyMixedMBR(CSTR): + _N_ins = 1 + _N_outs = 2 # [0] filtrate, [1] pumped flow + _outs_size_is_fixed = False + + def __init__(self, ID='', ins=None, outs=(), thermo=None, + init_with='WasteStream', isdynamic=True, + pumped_flow=50, solids_capture_rate=0.999, + V_max=1000, crossflow_air=None, + **kwargs): + super().__init__(ID, ins, outs, split=None, thermo=thermo, + init_with=init_with, V_max=V_max, isdynamic=isdynamic, + **kwargs) + self.pumped_flow = pumped_flow + self.solids_capture_rate = solids_capture_rate + self.crossflow_air = crossflow_air + + @property + def pumped_flow(self): + '''[float] Pumped flow rate, in m3/d''' + return self._Q_pump + @pumped_flow.setter + def pumped_flow(self, Q): + self._Q_pump = Q + + @property + def solids_capture_rate(self): + '''[float] Membrane solid capture rate, i.e., + filtrate-to-internal solids concentration ratio, unitless.''' + return self._f_rtn + @solids_capture_rate.setter + def solids_capture_rate(self, f): + if f < 0 or f > 1: + raise ValueError(f'membrane solids capture rate must be within [0,1], not {f}') + self._f_rtn = f + cmps = self._mixed.components + self._flt2in_conc_ratio = (1-cmps.x) + cmps.x * (1-f) + + @property + def crossflow_air(self): + '''[:class:`qsdsan.Process` or NoneType] + Membrane cross flow air specified for process modeling, such as `qsdsan.processes.DiffusedAeration`. + Ignored if DO setpoint is specified by the `aeration` attribute. + ''' + return self._cfa + @crossflow_air.setter + def crossflow_air(self, cfa): + if cfa is None or isinstance(cfa, Process): + self._cfa = cfa + else: + raise TypeError('crossflow_air must be a `Process` object or None, ' + f'not {type(cfa)}') + + split = property(CSTR.split.fget) + split.fset = None + + def _run(self): + '''Only to converge volumetric flows.''' + mixed = self._mixed + mixed.mix_from(self.ins) + cmps = mixed.components + Q = mixed.F_vol*24 # m3/d + Qp = self._Q_pump + f_rtn = self._f_rtn + xsplit = Qp / ((1-f_rtn)*(Q-Qp) + Qp) # mass split of solids to pumped flow + qsplit = Qp / Q + flt, rtn = self.outs + mixed.split_to(rtn, flt, xsplit*cmps.x + qsplit*(1-cmps.x)) + + def _compile_ODE(self): + aer = self._aeration + cfa = self._cfa + isa = isinstance + cmps = self.components + if self._model is None: + warn(f'{self.ID} was initialized without a suspended growth model, ' + f'and thus run as a non-reactive unit') + r = lambda state_arr: np.zeros(cmps.size) + else: + r = self._model.production_rates_eval + + _dstate = self._dstate + _update_dstate = self._update_dstate + V = self._V_max + Qp = self.pumped_flow + f_rtn = self.solids_capture_rate + xarr = cmps.x + gstrip = self.gas_stripping + if gstrip: + gas_idx = self.components.indices(self.gas_IDs) + if isa(aer, Process): kLa = aer.kLa + else: kLa = 0. + if cfa: kLa += cfa.kLa + S_gas_air = np.asarray(self.K_Henry)*np.asarray(self.p_gas_atm) + kLa_stripping = np.maximum(kLa*self.D_gas/self._D_O2, self.stripping_kLa_min) + hasexo = bool(len(self._exovars)) + f_exovars = self.eval_exo_dynamic_vars + + if isa(aer, (float, int)): + i = cmps.index(self._DO_ID) + def dy_dt(t, QC_ins, QC, dQC_ins): + QC[i] = aer + dydt_mbr(QC_ins, QC, V, Qp, f_rtn, xarr, _dstate) + if hasexo: QC = np.append(QC, f_exovars(t)) + _dstate[:-1] += r(QC) + if gstrip: _dstate[gas_idx] -= kLa_stripping * (QC[gas_idx] - S_gas_air) + _dstate[i] = 0 + _update_dstate() + else: + if cfa: + cfa_stoi = cfa._stoichiometry + cfa_frho = cfa.rate_function + dy_cfa = lambda QC: cfa_stoi * cfa_frho(QC) + else: + dy_cfa = lambda QC: 0. + + if isa(aer, Process): + aer_stoi = aer._stoichiometry + aer_frho = aer.rate_function + dy_aer = lambda QC: aer_stoi * aer_frho(QC) + else: + dy_aer = lambda QC: 0. + + def dy_dt(t, QC_ins, QC, dQC_ins): + dydt_mbr(QC_ins, QC, V, Qp, f_rtn, xarr, _dstate) + if hasexo: QC = np.append(QC, f_exovars(t)) + _dstate[:-1] += r(QC) + dy_aer(QC) + dy_cfa(QC) + if gstrip: _dstate[gas_idx] -= kLa_stripping * (QC[gas_idx] - S_gas_air) + _update_dstate() + self._ODE = dy_dt + + def _update_state(self): + arr = self._state + arr[arr < 1e-16] = 0. + arr[-1] = sum(ws.state[-1] for ws in self.ins) + flt, rtn = self.outs + Qp = self.pumped_flow + flt.state[:-1] = arr[:-1] * self._flt2in_conc_ratio + flt.state[-1] = arr[-1] - Qp + rtn.state[:-1] = arr[:-1] + rtn.state[-1] = Qp + + def _update_dstate(self): + arr = self._dstate + arr[-1] = sum(ws.dstate[-1] for ws in self.ins) + flt, rtn = self.outs + flt.dstate[:-1] = arr[:-1] * self._flt2in_conc_ratio + flt.dstate[-1] = arr[-1] + rtn.dstate[:-1] = arr[:-1] + rtn.dstate[-1] = 0 + + +#%% + +class PlugFlowMBR(PFR): + pass \ No newline at end of file From a37da6b8c0ff62dbf5859372512b4c97a5657e72 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 19 Nov 2024 14:15:58 -0800 Subject: [PATCH 04/14] add `CompletelyMixedMBR` doctest --- qsdsan/sanunits/_membrane_bioreactor.py | 85 ++++++++++++++++++++++++- 1 file changed, 82 insertions(+), 3 deletions(-) diff --git a/qsdsan/sanunits/_membrane_bioreactor.py b/qsdsan/sanunits/_membrane_bioreactor.py index 65e07212..43d1ff6d 100644 --- a/qsdsan/sanunits/_membrane_bioreactor.py +++ b/qsdsan/sanunits/_membrane_bioreactor.py @@ -1362,9 +1362,85 @@ def dydt_mbr(QC_ins, QC, V, Qp, f_rtn, xarr, _dstate): _dstate[:-1] = (Q_ins @ C_ins - (Q*(1-xarr) + (Qp+(Q-Qp)*(1-f_rtn))*xarr)*C)/V class CompletelyMixedMBR(CSTR): + ''' + Completely mixed membrane bioreactor, equivalent to a CSTR with ideal + membrane filtration at the outlet. + + See Also + -------- + :class:`qsdsan.processes.DiffusedAeration` + :class:`qsdsan.sanunits.CSTR` + + Examples + -------- + >>> from qsdsan import WasteStream, processes as pc, sanunits as su + >>> cmps = pc.create_asm1_cmps() + >>> ws = WasteStream('ws', H2O=100000, X_I=5, X_S=2, S_I=6) + >>> M1 = su.CompletelyMixedMBR('M1', ins=ws, pumped_flow=50, + ... solids_capture_rate=0.999, + ... V_max=1000, DO_ID='S_O') + >>> M1.simulate(t_span=(0,100), method='BDF') + >>> M1.show() + CompletelyMixedMBR: M1 + ins... + [0] ws + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (g/hr): S_I 6e+03 + X_I 5e+03 + X_S 2e+03 + H2O 1e+08 + WasteStream-specific properties: + pH : 7.0 + Alkalinity : 2.5 mg/L + COD : 129.6 mg/L + BOD : 11.6 mg/L + TC : 41.5 mg/L + TOC : 41.5 mg/L + TN : 3.0 mg/L + TP : 1.3 mg/L + TSS : 38.8 mg/L + outs... + [0] ws1 + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (g/hr): S_I 5.88e+03 + X_I 224 + X_S 89.6 + S_O 196 + H2O 9.79e+07 + WasteStream-specific properties: + pH : 7.0 + COD : 63.0 mg/L + BOD : 0.5 mg/L + TC : 20.2 mg/L + TOC : 20.2 mg/L + TN : 0.1 mg/L + TP : 0.6 mg/L + TSS : 1.8 mg/L + [1] ws2 + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (g/hr): S_I 125 + X_I 4.75e+03 + X_S 1.9e+03 + S_O 4.17 + H2O 2.07e+06 + WasteStream-specific properties: + pH : 7.0 + COD : 3253.1 mg/L + BOD : 529.2 mg/L + TC : 1041.0 mg/L + TOC : 1041.0 mg/L + TN : 136.9 mg/L + TP : 32.5 mg/L + TSS : 1774.0 mg/L + + >>> flt, rtn = M1.outs + >>> flt.get_TSS() / rtn.get_TSS() + 0.001 + + ''' _N_ins = 1 _N_outs = 2 # [0] filtrate, [1] pumped flow - _outs_size_is_fixed = False + _outs_size_is_fixed = True def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', isdynamic=True, @@ -1414,8 +1490,7 @@ def crossflow_air(self, cfa): raise TypeError('crossflow_air must be a `Process` object or None, ' f'not {type(cfa)}') - split = property(CSTR.split.fget) - split.fset = None + split = None def _run(self): '''Only to converge volumetric flows.''' @@ -1496,6 +1571,10 @@ def _update_state(self): arr = self._state arr[arr < 1e-16] = 0. arr[-1] = sum(ws.state[-1] for ws in self.ins) + for ws in self.outs: + if ws.state is None: + ws.state = np.zeros_like(arr) + ws.dstate = np.zeros_like(arr) flt, rtn = self.outs Qp = self.pumped_flow flt.state[:-1] = arr[:-1] * self._flt2in_conc_ratio From 469fe170bfe79e17e63660df3ee288902e5feaa1 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 19 Nov 2024 14:57:13 -0800 Subject: [PATCH 05/14] update `CompletelyMixedMBR` doctest --- qsdsan/sanunits/_membrane_bioreactor.py | 127 +++++++++++++++--------- 1 file changed, 82 insertions(+), 45 deletions(-) diff --git a/qsdsan/sanunits/_membrane_bioreactor.py b/qsdsan/sanunits/_membrane_bioreactor.py index 43d1ff6d..8a0149ac 100644 --- a/qsdsan/sanunits/_membrane_bioreactor.py +++ b/qsdsan/sanunits/_membrane_bioreactor.py @@ -1375,67 +1375,104 @@ class CompletelyMixedMBR(CSTR): -------- >>> from qsdsan import WasteStream, processes as pc, sanunits as su >>> cmps = pc.create_asm1_cmps() - >>> ws = WasteStream('ws', H2O=100000, X_I=5, X_S=2, S_I=6) - >>> M1 = su.CompletelyMixedMBR('M1', ins=ws, pumped_flow=50, - ... solids_capture_rate=0.999, - ... V_max=1000, DO_ID='S_O') - >>> M1.simulate(t_span=(0,100), method='BDF') + >>> ww = WasteStream('ww') + >>> ww.set_flow_by_concentration( + ... flow_tot=2000, + ... concentrations=dict(S_I=21.6, S_S=86.4, X_I=32.4, X_S=129.6, + ... S_NH=25, S_ND=2.78, X_ND=6.28, S_ALK=84), + ... units=('m3/d', 'mg/L')) + >>> asm = pc.ASM1() + >>> M1 = su.CompletelyMixedMBR('M1', ins=ww, outs=['filtrate', 'pumped'], + ... V_max=400, pumped_flow=50, solids_capture_rate=0.999, + ... aeration=2.0, DO_ID='S_O', + ... suspended_growth_model=asm) + >>> M1.set_init_conc(X_I=1000, S_I=30, S_S=5.0, X_S=100, X_BH=500, + ... X_BA=100, X_P=100, S_O=2, S_NH=2, S_ND=1, + ... X_ND=1, S_NO=20, S_ALK=84) + >>> M1.simulate(t_span=(0,400), method='BDF') >>> M1.show() CompletelyMixedMBR: M1 ins... - [0] ws + [0] ww phase: 'l', T: 298.15 K, P: 101325 Pa - flow (g/hr): S_I 6e+03 - X_I 5e+03 - X_S 2e+03 - H2O 1e+08 + flow (g/hr): S_I 1.8e+03 + S_S 7.2e+03 + X_I 2.7e+03 + X_S 1.08e+04 + S_NH 2.08e+03 + S_ND 232 + X_ND 523 + S_ALK 7e+03 + H2O 8.31e+07 WasteStream-specific properties: pH : 7.0 Alkalinity : 2.5 mg/L - COD : 129.6 mg/L - BOD : 11.6 mg/L - TC : 41.5 mg/L - TOC : 41.5 mg/L - TN : 3.0 mg/L - TP : 1.3 mg/L - TSS : 38.8 mg/L + COD : 270.0 mg/L + BOD : 137.1 mg/L + TC : 170.4 mg/L + TOC : 86.4 mg/L + TN : 36.0 mg/L + TP : 2.7 mg/L + TSS : 121.5 mg/L outs... - [0] ws1 + [0] filtrate phase: 'l', T: 298.15 K, P: 101325 Pa - flow (g/hr): S_I 5.88e+03 - X_I 224 - X_S 89.6 - S_O 196 - H2O 9.79e+07 + flow (g/hr): S_I 1.75e+03 + S_S 277 + X_I 101 + X_S 4.83 + X_BH 236 + X_BA 13.7 + X_P 44.1 + S_O 162 + S_NO 1.8e+03 + S_NH 61.7 + S_ND 60 + X_ND 0.323 + S_ALK 3.6e+03 + S_N2 254 + H2O 8.1e+07 WasteStream-specific properties: pH : 7.0 - COD : 63.0 mg/L - BOD : 0.5 mg/L - TC : 20.2 mg/L - TOC : 20.2 mg/L - TN : 0.1 mg/L - TP : 0.6 mg/L - TSS : 1.8 mg/L - [1] ws2 + COD : 29.9 mg/L + BOD : 4.2 mg/L + TC : 54.0 mg/L + TOC : 9.7 mg/L + TN : 24.0 mg/L + TP : 0.3 mg/L + TK : 0.0 mg/L + TSS : 3.7 mg/L + [1] pumped phase: 'l', T: 298.15 K, P: 101325 Pa - flow (g/hr): S_I 125 - X_I 4.75e+03 - X_S 1.9e+03 - S_O 4.17 - H2O 2.07e+06 + flow (g/hr): S_I 45 + S_S 7.09 + X_I 2.6e+03 + X_S 124 + X_BH 6.06e+03 + X_BA 351 + X_P 1.13e+03 + S_O 4.17 + S_NO 46 + S_NH 1.58 + S_ND 1.54 + X_ND 8.29 + S_ALK 92.2 + S_N2 6.51 + H2O 2.07e+06 WasteStream-specific properties: pH : 7.0 - COD : 3253.1 mg/L - BOD : 529.2 mg/L - TC : 1041.0 mg/L - TOC : 1041.0 mg/L - TN : 136.9 mg/L - TP : 32.5 mg/L - TSS : 1774.0 mg/L + COD : 4950.0 mg/L + BOD : 1779.8 mg/L + TC : 1795.1 mg/L + TOC : 1750.8 mg/L + TN : 381.0 mg/L + TP : 77.2 mg/L + TK : 17.3 mg/L + TSS : 3693.8 mg/L >>> flt, rtn = M1.outs - >>> flt.get_TSS() / rtn.get_TSS() - 0.001 + >>> flt.get_TSS() / rtn.get_TSS() # doctest: +ELLIPSIS + 0.001... ''' _N_ins = 1 From b120586f5e00b551cac94bfaf31e6a44fd223e93 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 20 Nov 2024 12:52:05 -0800 Subject: [PATCH 06/14] process add-on for aerobic digestion --- qsdsan/processes/__init__.py | 5 +- qsdsan/processes/_aerobic_digestion_addon.py | 106 ++++++++++++++++++ qsdsan/sanunits/_membrane_bioreactor.py | 3 +- .../sanunits/_suspended_growth_bioreactor.py | 69 ++++++++++-- 4 files changed, 173 insertions(+), 10 deletions(-) create mode 100644 qsdsan/processes/_aerobic_digestion_addon.py diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index ba9beecb..65e861c8 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -67,6 +67,7 @@ def __init__(self): from ._decay import * from ._kinetic_reaction import * from ._pm2 import * +from ._aerobic_digestion_addon import * from . import ( _aeration, @@ -77,7 +78,8 @@ def __init__(self): # _madm1, _decay, _kinetic_reaction, - _pm2 + _pm2, + _aerobic_digestion_addon, ) __all__ = ( @@ -90,4 +92,5 @@ def __init__(self): *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, + *_aerobic_digestion_addon.__all__, ) diff --git a/qsdsan/processes/_aerobic_digestion_addon.py b/qsdsan/processes/_aerobic_digestion_addon.py new file mode 100644 index 00000000..009c70a6 --- /dev/null +++ b/qsdsan/processes/_aerobic_digestion_addon.py @@ -0,0 +1,106 @@ +# -*- 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 qsdsan import Process +from thermosteam import settings +_load_components = settings.get_default_chemicals + +__all__ = ('ASM_AeDigAddOn',) + +class ASM_AeDigAddOn(Process): + ''' + Creates a `Process` object representing the degradation of particulate + inert organic materials that typically occur in an aerobic digester. + Stoichiometry is determined by rules of element conservation in corresponding + activated sludge models. + + Parameters + ---------- + k_dig : float, optional + The 1st-order degradation rate constant, in d^(-1). The default is 0.04. + + See Also + -------- + :class:`qsdsan.processes.ASM1` + :class:`qsdsan.processes.ASM2d` + :class:`qsdsan.processes.mASM2d` + + Examples + -------- + >>> import qsdsan.processes as pc + >>> cmps_asm1 = pc.create_asm1_cmps() + >>> dig_asm1 = pc.ASM_AeDigAddOn('dig_asm1') + >>> dig_asm1.show() + Process: dig_asm1 + [stoichiometry] X_I: -1.00 + X_S: 1.00 + S_NH: 0.0600 + S_ALK: 0.0515 + [reference] X_I + [rate equation] X_I*k_dig + [parameters] k_dig: 0.04 + [dynamic parameters] + + >>> cmps_masm2d = pc.create_masm2d_cmps(set_thermo=False) + >>> dig_masm2d = pc.ASM_AeDigAddOn('dig_masm2d', components=cmps_masm2d) + >>> dig_masm2d.show() + Process: dig_masm2d + [stoichiometry] S_NH4: 0.0265 + S_PO4: 0.000900 + S_IC: 0.0434 + X_I: -1.00 + X_S: 1.00 + [reference] X_I + [rate equation] X_I*k_dig + [parameters] k_dig: 0.04 + [dynamic parameters] + + ''' + + def __init__(self, ID, k_dig=0.04, components=None): + cmps = _load_components(components) + rxn = 'X_I -> X_S' + consrv = [] + if 'S_ALK' in cmps.IDs: + consrv.append('charge') + rxn += ' + [?]S_ALK' + elif 'S_IC' in cmps.IDs: + consrv.append('C') + rxn += ' + [?]S_IC' + + if 'S_NH' in cmps.IDs: + consrv.append('N') + rxn += ' + [?]S_NH' + elif 'S_NH4' in cmps.IDs: + consrv.append('N') + rxn += ' + [?]S_NH4' + + if 'S_PO4' in cmps.IDs: + consrv.append('P') + rxn += ' +[?]S_PO4' + + super().__init__(ID=ID, reaction=rxn, + rate_equation='k_dig*X_I', + ref_component='X_I', + components=cmps, + conserved_for=consrv, + parameters=('k_dig',)) + self.k_dig=k_dig + + @property + def k_dig(self): + '''[float] Degradation rate constant, in d^(-1).''' + return self._k + @k_dig.setter + def k_dig(self, k): + self._k = k + self.set_parameters(k_dig=k) \ No newline at end of file diff --git a/qsdsan/sanunits/_membrane_bioreactor.py b/qsdsan/sanunits/_membrane_bioreactor.py index 8a0149ac..cd401bab 100644 --- a/qsdsan/sanunits/_membrane_bioreactor.py +++ b/qsdsan/sanunits/_membrane_bioreactor.py @@ -35,7 +35,8 @@ __all__ = ('AnMBR', 'CompletelyMixedMBR', - 'PlugFlowMBR',) + # 'PlugFlowMBR', + ) #%% degassing = SanStream.degassing diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 4b8cc3e8..10a045e2 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -24,6 +24,7 @@ 'BatchExperiment', # 'SBR', 'PFR', + 'AerobicDigester', ) # def _add_aeration_to_growth_model(aer, model): @@ -332,19 +333,20 @@ def ODE(self): self._compile_ODE() return self._ODE - def _compile_ODE(self): - isa = isinstance - cmps = self.components - m = cmps.size - aer = self._aeration + def _init_model(self): if self._model is None: warn(f'{self.ID} was initialized without a suspended growth model, ' f'and thus run as a non-reactive unit') - r = lambda state_arr: np.zeros(m) - + r = lambda state_arr: np.zeros(self.components.size) else: # processes = _add_aeration_to_growth_model(aer, self._model) r = self._model.production_rates_eval + return r + + def _compile_ODE(self): + isa = isinstance + aer = self._aeration + r = self._init_model() _dstate = self._dstate _update_dstate = self._update_dstate @@ -1286,4 +1288,55 @@ def get_retained_mass(self, biomass_IDs): return mass @ self.V_tanks def _design(self): - pass \ No newline at end of file + pass + +#%% +from ..processes import ASM_AeDigAddOn + +class AerobicDigester(CSTR): + + def __init__(self, ID='', ins=None, outs=(), thermo=None, + init_with='WasteStream', V_max=1000, activated_sludge_model=None, + organic_particulate_inert_degradation_process=None, + aeration=1.0, DO_ID='S_O2', isdynamic=True, **kwargs): + super().__init__(ID, ins, outs, thermo=thermo, init_with=init_with, + V_max=V_max, aeration=aeration, DO_ID=DO_ID, + suspended_growth_model=activated_sludge_model, + isdynamic=isdynamic, **kwargs) + self.organic_particulate_inert_degradation_process = organic_particulate_inert_degradation_process + + @property + def organic_particulate_inert_degradation_process(self): + '''[:class:`Process` or NoneType] Process object for degradation of + particulate inert organic materials in the aerobic digester. If none + specified, will attempt to create a Process model according to components + by default.''' + return self._dig_addon + @organic_particulate_inert_degradation_process.setter + def organic_particulate_inert_degradation_process(self, proc): + if isinstance(proc, Process): + self._dig_addon = proc + elif proc is None: + if self._model is None: self._dig_addon = None + else: + ID = self._model.ID + '_particulate_inert_degrade' + self._dig_addon = ASM_AeDigAddOn( + ID=ID, + components=self.thermo.chemicals + ) + else: + raise TypeError('organic_particulate_inert_degradation_process must be' + f' a `Process` object if not None, not {type(proc)}') + + def _init_model(self): + if self._model is None: + warn(f'{self.ID} was initialized without an activated sludge model, ' + f'and thus run as a non-reactive unit') + r = lambda state_arr: np.zeros(self.components.size) + else: + dig = self.organic_particulate_inert_degradation_process + dig_stoi = dig._stoichiometry + dig_frho = dig.rate_function + asm_frate = self._model.production_rates_eval + r = lambda state_arr: asm_frate(state_arr) + dig_stoi * dig_frho(state_arr) + return r \ No newline at end of file From 0310be2d6fac03592ec6273f811bb432524ce4f2 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Fri, 22 Nov 2024 12:29:18 -0800 Subject: [PATCH 07/14] added doctest for `AerobicDigester` --- qsdsan/processes/_aerobic_digestion_addon.py | 3 +- .../sanunits/_suspended_growth_bioreactor.py | 99 ++++++++++++++++++- 2 files changed, 100 insertions(+), 2 deletions(-) diff --git a/qsdsan/processes/_aerobic_digestion_addon.py b/qsdsan/processes/_aerobic_digestion_addon.py index 009c70a6..440b604e 100644 --- a/qsdsan/processes/_aerobic_digestion_addon.py +++ b/qsdsan/processes/_aerobic_digestion_addon.py @@ -9,7 +9,7 @@ Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt for license details. ''' - +import numpy as np from qsdsan import Process from thermosteam import settings _load_components = settings.get_default_chemicals @@ -95,6 +95,7 @@ def __init__(self, ID, k_dig=0.04, components=None): conserved_for=consrv, parameters=('k_dig',)) self.k_dig=k_dig + self._stoichiometry = np.asarray(self._stoichiometry, dtype=float) @property def k_dig(self): diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 10a045e2..bc567109 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -1294,7 +1294,104 @@ def _design(self): from ..processes import ASM_AeDigAddOn class AerobicDigester(CSTR): + """ + Models an aerobic digester with activated sludge models, is a subclass of CSTR. + An additional degradation process of particulate inert organic materials is + considered in addition to typical activated sludge processes. + + Parameters + ---------- + activated_sludge_model : :class:`CompiledProcesses`, optional + The activated sludge model used for the biochemical reactions. The default is None. + organic_particulate_inert_degradation_process : :class:`Process`, optional + The degradation process of the particulate inert organic materials. The default is None. + See Also + -------- + :class:`qsdsan.processes.ASM1` + :class:`qsdsan.processes.ASM2d` + :class:`qsdsan.processes.mASM2d` + :class:`qsdsan.processes.ASM_AeDigAddOn` + :class:`qsdsan.sanunites.CSTR` + + Examples + -------- + >>> from qsdsan import WasteStream, processes as pc, sanunits as su + >>> cmps = pc.create_asm1_cmps() + >>> twas = WasteStream('thickened_WAS') + >>> twas.set_flow_by_concentration( + ... flow_tot=50, + ... concentrations=dict( + ... S_I=30, S_S=1, X_I=17000, X_S=800, X_BH=38000, X_BA=2300, + ... X_P=6500, S_O=0.5, S_NH=2, S_ND=1, X_ND=65, S_NO=10, + ... S_N2=25, S_ALK=84), + ... units=('m3/d', 'mg/L')) + >>> asm = pc.ASM1() + >>> AED = su.AerobicDigester('AED', ins=twas, outs=('digested_WAS',), + ... V_max=3000, activated_sludge_model=asm, + ... DO_ID='S_O', aeration=1.0) + >>> AED.simulate(t_span=(0, 400), method='BDF') + >>> AED.show() + AerobicDigester: AED + ins... + [0] thickened_WAS + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (g/hr): S_I 62.5 + S_S 2.08 + X_I 3.54e+04 + X_S 1.67e+03 + X_BH 7.92e+04 + X_BA 4.79e+03 + X_P 1.35e+04 + S_O 1.04 + S_NO 20.8 + S_NH 4.17 + S_ND 2.08 + X_ND 135 + S_ALK 175 + S_N2 52.1 + H2O 1.99e+06 + WasteStream-specific properties: + pH : 7.0 + Alkalinity : 2.5 mg/L + COD : 64631.0 mg/L + BOD : 23300.3 mg/L + TC : 22923.3 mg/L + TOC : 22839.3 mg/L + TN : 4712.0 mg/L + TP : 1009.0 mg/L + TK : 223.2 mg/L + TSS : 48450.0 mg/L + outs... + [0] digested_WAS + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (g/hr): S_I 62.5 + S_S 3.58e+04 + X_I 1.04e+04 + X_S 123 + X_BH 9.6e+03 + X_BA 1.59e+03 + X_P 2.77e+04 + S_O 2.08 + S_NO 4.17e+03 + S_NH 0.101 + S_ND 0.975 + X_ND 8.74 + S_N2 2.51e+03 + H2O 2e+06 + WasteStream-specific properties: + pH : 7.0 + Alkalinity : 2.5 mg/L + COD : 40987.6 mg/L + BOD : 15416.1 mg/L + TC : 13985.1 mg/L + TOC : 13985.1 mg/L + TN : 3534.1 mg/L + TP : 458.2 mg/L + TK : 89.3 mg/L + TSS : 17813.2 mg/L + + """ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', V_max=1000, activated_sludge_model=None, organic_particulate_inert_degradation_process=None, @@ -1319,7 +1416,7 @@ def organic_particulate_inert_degradation_process(self, proc): elif proc is None: if self._model is None: self._dig_addon = None else: - ID = self._model.ID + '_particulate_inert_degrade' + ID = self._model.__class__.__name__ + '_particulate_inert_degrade' self._dig_addon = ASM_AeDigAddOn( ID=ID, components=self.thermo.chemicals From 69828f73d9068d4563928058378e5bb1efe1baf2 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Fri, 22 Nov 2024 12:55:13 -0800 Subject: [PATCH 08/14] Update _aerobic_digestion_addon.py --- qsdsan/processes/_aerobic_digestion_addon.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/qsdsan/processes/_aerobic_digestion_addon.py b/qsdsan/processes/_aerobic_digestion_addon.py index 440b604e..bfdd267a 100644 --- a/qsdsan/processes/_aerobic_digestion_addon.py +++ b/qsdsan/processes/_aerobic_digestion_addon.py @@ -41,10 +41,10 @@ class ASM_AeDigAddOn(Process): >>> dig_asm1 = pc.ASM_AeDigAddOn('dig_asm1') >>> dig_asm1.show() Process: dig_asm1 - [stoichiometry] X_I: -1.00 - X_S: 1.00 - S_NH: 0.0600 - S_ALK: 0.0515 + [stoichiometry] X_I: -1 + X_S: 1 + S_NH: 0.06 + S_ALK: 0.0514 [reference] X_I [rate equation] X_I*k_dig [parameters] k_dig: 0.04 @@ -55,14 +55,14 @@ class ASM_AeDigAddOn(Process): >>> dig_masm2d.show() Process: dig_masm2d [stoichiometry] S_NH4: 0.0265 - S_PO4: 0.000900 - S_IC: 0.0434 - X_I: -1.00 - X_S: 1.00 + S_PO4: 0.0009 + S_IC: 0.0433 + X_I: -1 + X_S: 1 [reference] X_I [rate equation] X_I*k_dig [parameters] k_dig: 0.04 - [dynamic parameters] + [dynamic parameters] ''' From 82524ac09716cb24d8fbba248f7c81e3d48f3945 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Mon, 25 Nov 2024 13:22:05 -0800 Subject: [PATCH 09/14] A simple influent model for mASM2d --- qsdsan/processes/_asm2d.py | 115 ++++++++++++++++++++++++++++++++++++- 1 file changed, 113 insertions(+), 2 deletions(-) diff --git a/qsdsan/processes/_asm2d.py b/qsdsan/processes/_asm2d.py index c2a99b7e..3783b80e 100644 --- a/qsdsan/processes/_asm2d.py +++ b/qsdsan/processes/_asm2d.py @@ -12,7 +12,7 @@ import numpy as np from thermosteam.utils import chemicals_user from thermosteam import settings -from qsdsan import Component, Components, Processes, CompiledProcesses +from qsdsan import Component, Components, WasteStream, Processes, CompiledProcesses from ..utils import ospath, data_path, load_data from . import Monod, ion_speciation from scipy.optimize import brenth @@ -20,7 +20,8 @@ __all__ = ('create_asm2d_cmps', 'ASM2d', - 'create_masm2d_cmps', 'mASM2d') + 'create_masm2d_cmps', 'mASM2d', + 'create_masm2d_inf') _path = ospath.join(data_path, 'process_data/_asm2d.tsv') _load_components = settings.get_default_chemicals @@ -990,3 +991,113 @@ def set_pKsps(self, ps): K *= m2m**abs(v) Ksp_mass.append(K) self.rate_function._params['Ksp'] = np.array(Ksp_mass) + +#%% + +def create_masm2d_inf( + ID, Q, Q_unit='m3/d', T=298.15, + COD=430, NH4_N=25.0, PO4_P=8.0, alkalinity=7.0, + fr_SI=0.05, fr_SF=0.2, fr_SA=0.0, + fr_XI=0.13, fr_XH=0.0, fr_XAUT=0.0, + fr_XPAO=0.0, fr_XPHA=0.0, X_PP=0.0, + S_NO3=0.0, S_O2=0.0, S_N2=18, + S_Ca=140, S_Mg=50, S_K=28, S_Na=87, S_Cl=425, + X_AlOH=0, X_AlPO4=0, X_FeOH=0, X_FePO4=0, + X_CaCO3=0, X_ACP=0, X_MgCO3=0, X_newb=0, X_struv=0 + ): + ''' + Convenient function to create an influent `WasteStream` object with `mASM2d` + state variables based on specified bulk properties. + + Parameters + ---------- + ID : str + Unique identification for the `WasteStream` object. + Q : float + Total volumetric flow rate. + Q_unit : str, optional + Unit of measurement for flow rate. The default is 'm3/d'. + T : float, optional + Temperature, in K. The default is 298.15. + COD : float, optional + Total chemical oxygen demand, not accounting for electron acceptors like + dissvoled O2 or nitrate/nitrite, in mg-COD/L. The default is 430. + NH4_N : float, optional + Ammonium nitrogen, in mg-N/L. The default is 25.0. + PO4_P : float, optional + Ortho-phosphate, in mg-P/L. The default is 8.0. + alkalinity : float, optional + In mmol/L. The default is 7.0. + fr_SI : float, optional + Soluble inert fraction of total COD. The default is 0.05. + fr_SF : float, optional + Fermentable biodegradable fraction of total COD. The default is 0.2. + fr_SA : float, optional + VFA fraction of total COD. The default is 0.0. + fr_XI : float, optional + Particulate inert fraction of total COD. The default is 0.13. + fr_XH : float, optional + Heterotrophic biomass fraction of total COD. The default is 0.0. + fr_XAUT : float, optional + Autotrophic biomass fraction of total COD. The default is 0.0. + fr_XPAO : float, optional + Phosphorus accumulating biomass fraction of total COD. The default is 0.0. + fr_XPHA : float, optional + PHA fraction of total COD. The default is 0.0. + X_PP : float, optional + Poly-phosphate in mg-P/L. The default is 0.0. + S_NO3 : float, optional + Nitrate and nitrite in mg-N/L. The default is 0.0. + S_O2 : float, optional + Dissolved oxygen in mg-O2/L. The default is 0.0. + S_N2 : float, optional + Dissolved nitrogen gas in mg-N/L. The default is 18. + S_Ca : float, optional + Total soluble calcium in mg-Ca/L. The default is 140. + S_Mg : float, optional + Total soluble magnesium in mg-Mg/L. The default is 50. + S_K : float, optional + Total soluble potassium in mg-K/L. The default is 28. + S_Na : float, optional + Other cation, in mg-Na/L. The default is 87. + S_Cl : float, optional + Other anion, in mg-Cl/L. The default is 425. + X_AlOH : float, optional + Aluminum hydroxide [mg/L]. The default is 0. + X_AlPO4 : float, optional + Aluminum phosphate [mg/L]. The default is 0. + X_FeOH : float, optional + Iron hydroxide [mg/L]. The default is 0. + X_FePO4 : float, optional + Iron phosphate [mg/L]. The default is 0. + X_CaCO3 : float, optional + Calcium carbonate [mg/L]. The default is 0. + X_ACP : float, optional + Calcium phosphate [mg/L]. The default is 0. + X_MgCO3 : float, optional + Magnesium carbonate [mg/L]. The default is 0. + X_newb : float, optional + Newbryite [mg/L]. The default is 0. + X_struv : float, optional + Struvite [mg/L]. The default is 0. + + ''' + + fr_xs = 1.0-fr_SI-fr_SF-fr_SA-fr_XI-fr_XH-fr_XAUT-fr_XPAO-fr_XPHA + if fr_xs < 0: + raise ValueError('The sum of all COD fractions of organic materials must ' + 'not exceed 1.') + + inf = WasteStream(ID, T=T, SAlk=alkalinity) + concs = dict( + S_NH4=NH4_N, S_PO4=PO4_P, S_IC=alkalinity*12, + S_I=COD*fr_SI, S_F=COD*fr_SF, S_A=COD*fr_SA, + X_S=COD*fr_xs, X_I=COD*fr_XI, X_H=COD*fr_XH, + X_AUT=COD*fr_XAUT, X_PAO=COD*fr_XPAO, X_PHA=COD*fr_XPHA, + X_PP=X_PP, S_NO3=S_NO3, S_O2=S_O2, S_N2=S_N2, + S_Ca=S_Ca, S_Mg=S_Mg, S_K=S_K, S_Na=S_Na, S_Cl=S_Cl, + X_AlOH=X_AlOH, X_AlPO4=X_AlPO4, X_FeOH=X_FeOH, X_FePO4=X_FePO4, + X_CaCO3=X_CaCO3, X_ACP=X_ACP, X_MgCO3=X_MgCO3, X_newb=X_newb, X_struv=X_struv + ) + inf.set_flow_by_concentration(Q, concs, (Q_unit, 'mg/L')) + return inf From b714a90d28db44854622a1d4bee85f42386f3242 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Mon, 25 Nov 2024 13:46:46 -0800 Subject: [PATCH 10/14] minor bug fix --- qsdsan/_waste_stream.py | 10 +++++----- qsdsan/sanunits/_clarifier.py | 4 ++-- qsdsan/sanunits/_electrochemical_cell.py | 2 +- qsdsan/sanunits/_membrane_bioreactor.py | 2 +- qsdsan/sanunits/_sludge_treatment.py | 4 ++-- qsdsan/sanunits/_suspended_growth_bioreactor.py | 6 +++--- 6 files changed, 14 insertions(+), 14 deletions(-) diff --git a/qsdsan/_waste_stream.py b/qsdsan/_waste_stream.py index ca737b98..37c5843d 100644 --- a/qsdsan/_waste_stream.py +++ b/qsdsan/_waste_stream.py @@ -450,7 +450,7 @@ def _wastestream_info(self, details=True, concentrations=None, N=15): _ws_info += '\n' # Only non-zero properties are shown _ws_info += int(bool(self.pH))*f' pH : {self.pH:.1f}\n' - _ws_info += int(bool(self.SAlk))*f' Alkalinity : {self.SAlk:.1f} mg/L\n' + _ws_info += int(bool(self.SAlk))*f' Alkalinity : {self.SAlk:.1f} mmol/L\n' if details: _ws_info += int(bool(self.COD)) *f' COD : {self.COD:.1f} mg/L\n' _ws_info += int(bool(self.BOD)) *f' BOD : {self.BOD:.1f} mg/L\n' @@ -1299,7 +1299,7 @@ def codstates_inf_model(cls, ID='', flow_tot=0., units = ('L/hr', 'mg/L'), ... 9.96e+04 WasteStream-specific properties: pH : 7.0 - Alkalinity : 10.0 mg/L + Alkalinity : 10.0 mmol/L COD : 430.0 mg/L BOD : 221.8 mg/L TC : 265.0 mg/L @@ -1522,7 +1522,7 @@ def codbased_inf_model(cls, ID='', flow_tot=0., units = ('L/hr', 'mg/L'), ... 9.96e+04 WasteStream-specific properties: pH : 7.0 - Alkalinity : 10.0 mg/L + Alkalinity : 10.0 mmol/L COD : 430.0 mg/L BOD : 249.4 mg/L TC : 265.0 mg/L @@ -1751,7 +1751,7 @@ def bodbased_inf_model(cls, ID='', flow_tot=0., units = ('L/hr', 'mg/L'), ... 9.96e+04 WasteStream-specific properties: pH : 7.0 - Alkalinity : 10.0 mg/L + Alkalinity : 10.0 mmol/L COD : 431.0 mg/L BOD : 250.0 mg/L TC : 264.9 mg/L @@ -1983,7 +1983,7 @@ def sludge_inf_model(cls, ID='', flow_tot=0., units = ('L/hr', 'mg/L'), ... 9.88e+04 WasteStream-specific properties: pH : 7.0 - Alkalinity : 10.0 mg/L + Alkalinity : 10.0 mmol/L COD : 10814.4 mg/L BOD : 1744.3 mg/L TC : 4246.5 mg/L diff --git a/qsdsan/sanunits/_clarifier.py b/qsdsan/sanunits/_clarifier.py index b0053903..663dcce9 100644 --- a/qsdsan/sanunits/_clarifier.py +++ b/qsdsan/sanunits/_clarifier.py @@ -1049,7 +1049,7 @@ class PrimaryClarifierBSM2(SanUnit): H2O 1e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 23873.0 mg/L BOD : 14963.2 mg/L TC : 8298.3 mg/L @@ -1270,7 +1270,7 @@ class PrimaryClarifier(IdealClarifier): H2O 1e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 23873.0 mg/L BOD : 14963.2 mg/L TC : 8298.3 mg/L diff --git a/qsdsan/sanunits/_electrochemical_cell.py b/qsdsan/sanunits/_electrochemical_cell.py index 026d7046..d3c9fae3 100644 --- a/qsdsan/sanunits/_electrochemical_cell.py +++ b/qsdsan/sanunits/_electrochemical_cell.py @@ -95,7 +95,7 @@ class ElectrochemicalCell(SanUnit): O2 0.104 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L TC : 1860.0 mg/L TP : 31.9 mg/L TK : 1470.0 mg/L diff --git a/qsdsan/sanunits/_membrane_bioreactor.py b/qsdsan/sanunits/_membrane_bioreactor.py index cd401bab..92076b0c 100644 --- a/qsdsan/sanunits/_membrane_bioreactor.py +++ b/qsdsan/sanunits/_membrane_bioreactor.py @@ -1407,7 +1407,7 @@ class CompletelyMixedMBR(CSTR): H2O 8.31e+07 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 270.0 mg/L BOD : 137.1 mg/L TC : 170.4 mg/L diff --git a/qsdsan/sanunits/_sludge_treatment.py b/qsdsan/sanunits/_sludge_treatment.py index b075f2ed..c817830f 100644 --- a/qsdsan/sanunits/_sludge_treatment.py +++ b/qsdsan/sanunits/_sludge_treatment.py @@ -117,7 +117,7 @@ class Thickener(SanUnit): H2O 1e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 23873.0 mg/L BOD : 14963.2 mg/L TC : 8298.3 mg/L @@ -510,7 +510,7 @@ class Incinerator(SanUnit): H2O 1e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 23873.0 mg/L BOD : 14963.2 mg/L TC : 8298.3 mg/L diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index bc567109..6c5493a9 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -891,7 +891,7 @@ class PFR(SanUnit): H2O 1.53e+09 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 4389.1 mg/L BOD : 1563.3 mg/L TC : 1599.8 mg/L @@ -1353,7 +1353,7 @@ class AerobicDigester(CSTR): H2O 1.99e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 64631.0 mg/L BOD : 23300.3 mg/L TC : 22923.3 mg/L @@ -1381,7 +1381,7 @@ class AerobicDigester(CSTR): H2O 2e+06 WasteStream-specific properties: pH : 7.0 - Alkalinity : 2.5 mg/L + Alkalinity : 2.5 mmol/L COD : 40987.6 mg/L BOD : 15416.1 mg/L TC : 13985.1 mg/L From 38a791b4399973335c8cfb37995ab81e6e4ccd83 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 27 Nov 2024 14:56:56 -0800 Subject: [PATCH 11/14] minor update --- qsdsan/_process.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsdsan/_process.py b/qsdsan/_process.py index 5a5e768b..824e6ed1 100644 --- a/qsdsan/_process.py +++ b/qsdsan/_process.py @@ -1041,8 +1041,8 @@ def load_from_file(cls, path='', components=None, data=None, stoichio = proc[cmp_IDs] if data.columns[-1] in cmp_IDs: rate_eq = None else: - if pd.isna(proc[-1]): rate_eq = None - else: rate_eq = proc[-1] + if pd.isna(proc.iloc[-1]): rate_eq = None + else: rate_eq = proc.iloc[-1] stoichio = stoichio[-pd.isna(stoichio)].to_dict() ref = None for k,v in stoichio.items(): From 68a3ae1eee46dba51988c6ba20d043e20086f86c Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 3 Dec 2024 19:26:16 -0800 Subject: [PATCH 12/14] minor update to improve numerical stability --- qsdsan/sanunits/_junction.py | 4 ++++ qsdsan/sanunits/_suspended_growth_bioreactor.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/qsdsan/sanunits/_junction.py b/qsdsan/sanunits/_junction.py index f99809ed..60acce2c 100644 --- a/qsdsan/sanunits/_junction.py +++ b/qsdsan/sanunits/_junction.py @@ -2700,6 +2700,8 @@ def adm1p2masm2d(adm_vals): fraction_dissolve = max(0, min(1, - S_IC / xc_mmp)) asm_vals -= fraction_dissolve * X_CaCO3 * cac_sto asm_vals -= fraction_dissolve * X_MgCO3 * mgc_sto + if asm_vals[8] < 0: + asm_vals[8] = 0 if S_IN < 0: xn_mmp = sum(asm_vals[_mmp_idx] * mmp_in) if xn_mmp > 0: @@ -3022,6 +3024,8 @@ def masm2d2adm1p(asm_vals): fraction_dissolve = max(0, min(1, - S_IC / xc_mmp)) adm_vals -= fraction_dissolve * X_CaCO3 * cac_sto adm_vals -= fraction_dissolve * X_MgCO3 * mgc_sto + if adm_vals[9] < 0: + adm_vals[9] = 0 if S_IN < 0: xn_mmp = sum(adm_vals[_mmp_idx] * mmp_in) if xn_mmp > 0: diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 6c5493a9..19c6b2cf 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -1167,7 +1167,7 @@ def _init_state(self): def _update_state(self): out, = self.outs ncol = self._ncol - self._state[self._state < 2.2e-16] = 0. + self._state[self._state < 1e-16] = 0. self._state[self._Qs_idx] = self._Qs if out.state is None: out.state = np.zeros(ncol) out.state[:-1] = self._state[-ncol:-1] From c43e97948fc4bcfcb2a4c02425b514c4cf403a01 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 4 Dec 2024 12:40:48 -0800 Subject: [PATCH 13/14] fix adm1p algebraic H2 solver --- qsdsan/processes/_adm1_p_extension.py | 8 +++++--- qsdsan/sanunits/_junction.py | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py index aa4d1818..77f6305d 100644 --- a/qsdsan/processes/_adm1_p_extension.py +++ b/qsdsan/processes/_adm1_p_extension.py @@ -1037,7 +1037,7 @@ def __new__(cls, components=None, path=None, #!!! new parameter KS_IP*P_mw, np.array(k_mmp), Ksp_mass, np.array(K_dis), K_AlOH, K_FeOH])) - + 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[45] @@ -1045,6 +1045,7 @@ def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): stoichio = f_stoichio(state_arr) # should return the stoichiometric coefficients of S_h2 for all processes return Q/V_liq*(S_h2_in - S_h2) + np.dot(rxn, stoichio) + grad_rhosp = np.zeros(5) X_biop = np.zeros(5) def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): @@ -1060,12 +1061,13 @@ def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): X_biop[:] = state_arr[[18,19,19,20,22]] substrates = state_arr[2:6] - S_va, S_bu, S_IN = state_arr[[3,4,10]] + S_va, S_bu, S_IN, S_IP = state_arr[[3,4,10,11]] Iph = Hill_inhibit(h, pH_ULs, pH_LLs)[[2,3,4,5,7]] Iin = substr_inhibit(S_IN, KS_IN) + Iip = substr_inhibit(S_IP, KS_IP) grad_Ih2 = grad_non_compet_inhibit(S_h2, KIs_h2) - grad_rhosp[:] = ks * X_biop * Iph * Iin + grad_rhosp[:] = ks * X_biop * Iph * Iin * Iip grad_rhosp[:-1] *= substr_inhibit(substrates, Ks) * grad_Ih2 if S_va > 0: grad_rhosp[1] *= 1/(1+S_bu/S_va) if S_bu > 0: grad_rhosp[2] *= 1/(1+S_va/S_bu) diff --git a/qsdsan/sanunits/_junction.py b/qsdsan/sanunits/_junction.py index 60acce2c..56c4a0e8 100644 --- a/qsdsan/sanunits/_junction.py +++ b/qsdsan/sanunits/_junction.py @@ -2659,7 +2659,7 @@ def adm1p2masm2d(adm_vals): X_newb, X_ACP, X_MgCO3, X_AlOH, X_AlPO4, X_FeOH, X_FePO4, \ S_Na, S_Cl, H2O = _adm_vals - if S_h2 > 0 or S_ch4 > 0: warn('Ignored dissolved H2 or CH4.') + if S_h2 > 0 or S_ch4 > 0: warn('Ignored dissolved H2 or CH4 in ADM1p-to-mASM2d interface model.') S_NH4 = S_IN S_PO4 = S_IP From 2155dd8474647f94f2031f98811bbe444f7a64e4 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 4 Dec 2024 13:14:18 -0800 Subject: [PATCH 14/14] troubleshoot adm1p --- qsdsan/processes/_adm1_p_extension.py | 11 ++++++----- qsdsan/sanunits/_anaerobic_reactor.py | 7 ++++--- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py index 77f6305d..b850db4b 100644 --- a/qsdsan/processes/_adm1_p_extension.py +++ b/qsdsan/processes/_adm1_p_extension.py @@ -1038,7 +1038,7 @@ def __new__(cls, components=None, path=None, KS_IP*P_mw, np.array(k_mmp), Ksp_mass, np.array(K_dis), K_AlOH, K_FeOH])) - def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): + def adm1p_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): state_arr[7] = S_h2 Q = state_arr[45] rxn = _rhos_adm1p(state_arr, params, h=h) @@ -1048,7 +1048,7 @@ def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): grad_rhosp = np.zeros(5) X_biop = np.zeros(5) - def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): + def adm1p_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'][[5,6,7,8,10]] Ks = params['half_sat_coeffs'][2:6] @@ -1056,6 +1056,7 @@ def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): pH_ULs = params['pH_ULs'] pH_LLs = params['pH_LLs'] KS_IN = params['KS_IN'] + KS_IP = params['KS_IP'] KIs_h2 = params['KIs_h2'] kLa = params['kLa'] @@ -1076,12 +1077,12 @@ def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): stoichio = f_stoichio(state_arr) Q = state_arr[45] - return -Q/V_liq + np.dot(grad_rhos, stoichio[[5,6,7,8,10]]) + kLa*stoichio[-3] + return -Q/V_liq + np.dot(grad_rhosp, 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 + dct['dydt_Sh2_AD'] = adm1p_dydt_Sh2_AD + dct['grad_dydt_Sh2_AD'] = adm1p_grad_dydt_Sh2_AD return self def set_half_sat_K(self, K, process): diff --git a/qsdsan/sanunits/_anaerobic_reactor.py b/qsdsan/sanunits/_anaerobic_reactor.py index 496cb061..c800af47 100644 --- a/qsdsan/sanunits/_anaerobic_reactor.py +++ b/qsdsan/sanunits/_anaerobic_reactor.py @@ -583,8 +583,9 @@ def h2_stoichio(state_arr): dydt_Sh2_AD = self.model.dydt_Sh2_AD grad_dydt_Sh2_AD = self.model.grad_dydt_Sh2_AD def solve_h2(QC, S_in, T, h=h): - Ka = params['Ka_base'] * T_correction_factor(params['T_base'], T, params['Ka_dH']) - if h == None: h = solve_pH(QC, Ka, unit_conversion) + if h == None: + Ka = params['Ka_base'] * T_correction_factor(params['T_base'], T, params['Ka_dH']) + h = solve_pH(QC, Ka, unit_conversion) # S_h2_0 = QC[h2_idx] S_h2_0 = 2.8309E-07 S_h2_in = S_in[h2_idx] @@ -596,7 +597,7 @@ def solve_h2(QC, S_in, T, h=h): def update_h2_dstate(dstate): dstate[h2_idx] = 0. else: - solve_h2 = lambda QC, S_ins, T: QC[h2_idx] + solve_h2 = lambda QC, S_in, T: QC[h2_idx] def update_h2_dstate(dstate): pass def dy_dt(t, QC_ins, QC, dQC_ins):