From 2420d97390bc8e6e6071b275c98f3ee13cc27cc2 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Mon, 22 Apr 2024 17:44:23 -0500 Subject: [PATCH 01/22] added asm2d components to pm2 components finished working on rho --- qsdsan/data/process_data/_pm2_asm2d.tsv | 39 ++ qsdsan/processes/_pm2_asm2d.py | 793 ++++++++++++++++++++++++ 2 files changed, 832 insertions(+) create mode 100644 qsdsan/data/process_data/_pm2_asm2d.tsv create mode 100644 qsdsan/processes/_pm2_asm2d.py diff --git a/qsdsan/data/process_data/_pm2_asm2d.tsv b/qsdsan/data/process_data/_pm2_asm2d.tsv new file mode 100644 index 00000000..7bc0ff7f --- /dev/null +++ b/qsdsan/data/process_data/_pm2_asm2d.tsv @@ -0,0 +1,39 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +photoadaptation 1 +ammonium_uptake -1 1 +phosphorus_uptake -1 1 +growth_pho 1 ? 1 ? ? +carbohydrate_storage_pho 1 ? 1 +lipid_storage_pho 1 ? 1 +carbohydrate_growth_pho 1 (-Y_CH_PHO/Y_X_ALG_PHO) ? ? ? ? +lipid_growth_pho 1 (-Y_LI_PHO/Y_X_ALG_PHO) ? ? ? ? +carbohydrate_maintenance_pho -1 ? -1 +lipid_maintenance_pho -1 ? -1 +endogenous_respiration_pho -1 ? -1 ? ? +growth_ace 1 ? (-1)/Y_X_ALG_HET_ACE ? ? ? +carbohydrate_storage_ace 1 ? (-1)/Y_CH_ND_HET_ACE ? +lipid_storage_ace 1 ? (-1)/Y_LI_ND_HET_ACE ? +carbohydrate_growth_ace 1 (-Y_CH_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +lipid_growth_ace 1 (-Y_LI_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +carbohydrate_maintenance_ace -1 ? -1 +lipid_maintenance_ace -1 ? -1 +endogenous_respiration_ace -1 ? -1 ? ? +growth_glu 1 ? (-1)/Y_X_ALG_HET_GLU ? ? ? +carbohydrate_storage_glu 1 ? (-1)/Y_CH_ND_HET_GLU ? +lipid_storage_glu 1 ? (-1)/Y_LI_ND_HET_GLU ? +carbohydrate_growth_glu 1 (-Y_CH_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +lipid_growth_glu 1 (-Y_LI_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +carbohydrate_maintenance_glu -1 ? -1 +lipid_maintenance_glu -1 ? -1 +endogenous_respiration_glu -1 ? -1 ? ? +aero_hydrolysis 1-f_SI ? ? ? f_SI -1 +anox_hydrolysis 1-f_SI ? ? ? f_SI -1 +anae_hydrolysis 1-f_SI ? ? ? f_SI -1 +hetero_growth_S_F (-1)/Y_H 1-1/Y_H ? ? ? 1 +hetero_growth_S_A (-1)/Y_H 1-1/Y_H ? ? ? 1 +denitri_S_F (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +denitri_S_A (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +ferment 1 -1 ? ? ? +hetero_lysis ? ? ? f_XI_H 1-f_XI_H -1 +auto_aero_growth (Y_A-32/7)/Y_A ? 1/Y_A ? ? 1 +auto_lysis ? ? ? f_XI_AUT 1-f_XI_AUT -1 diff --git a/qsdsan/processes/_pm2_asm2d.py b/qsdsan/processes/_pm2_asm2d.py new file mode 100644 index 00000000..08a0a294 --- /dev/null +++ b/qsdsan/processes/_pm2_asm2d.py @@ -0,0 +1,793 @@ +# -*- coding: utf-8 -*- +''' +QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems + +This module is developed by: + Ga-Yeong Kim + Joy Zhang + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt +for license details. +''' +# from thermosteam.utils import chemicals_user +from thermosteam import settings +from qsdsan import Component, Components, Process, Processes, CompiledProcesses +from qsdsan.utils import ospath, data_path +import numpy as np + +__all__ = ('create_pm2_asm2d_cmps', 'PM2ASM2d') + +_path = ospath.join(data_path, 'process_data/_pm2_asm2d.tsv') +# _load_components = settings.get_default_chemicals + +#%% +# ============================================================================= +# PM2ASM2d-specific components +# ============================================================================= + +def create_pm2_asm2d_cmps(set_thermo=True): + cmps = Components.load_default() + + # X_CHL (g Chl/m^3) + X_CHL = Component(ID = 'X_CHL', + formula = 'C55H72MgN4O5', + description = 'Chlorophyll content of cells', + particle_size = 'Particulate', + degradability = 'Slowly', + organic = True) + + # X_ALG (g COD/m^3) + X_ALG = cmps.X_OHO.copy('X_ALG') + X_ALG.description = 'Concentration of carbon-accumulating mixotrophic organisms' + X_ALG.formula = 'CH1.8O0.5N0.2P0.018' + X_ALG.f_BOD5_COD = X_ALG.f_uBOD_COD = None + X_ALG.f_Vmass_Totmass = 0.89 + + # X_CH (g COD/m^3) + X_CH = cmps.X_GAO_Gly.copy('X_CH') + X_CH.description = 'Concentration of stored carbohydrates' + X_CH.formula = 'CH2O' + X_CH.f_BOD5_COD = X_CH.f_uBOD_COD = None + + # X_LI (g COD/m^3) + X_LI = cmps.X_GAO_Gly.copy('X_LI') + X_LI.description = 'Concentration of stored lipids' + X_LI.formula = 'CH1.92O0.118' + X_LI.f_BOD5_COD = X_LI.f_uBOD_COD = None + + # S_CO2 (g CO2/m^3) + S_CO2 = Component.from_chemical(ID = 'S_CO2', + chemical = 'CO2', + description = 'Soluble carbon dioxide', + particle_size = 'Soluble', + degradability = 'Undegradable', + organic = False) + + # S_A (g COD/m^3) + S_A = cmps.S_Ac.copy('S_A') + S_A.description = 'Concentration of extracellular dissolved organic carbon (acetate)' + + # S_F (g COD/m^3) + S_F = Component.from_chemical(ID = 'S_F', + chemical = 'glucose', + description = 'Concentration of extracellular dissolved organic carbon (glucose)', + measured_as = 'COD', + particle_size = 'Soluble', + degradability = 'Readily', + organic = True) + + # S_O2 (g O2/m^3) + S_O2 = cmps.S_O2.copy('S_O2') + S_O2.description = ('Concentration of dissolved oxygen') + + # S_NH (g N/m^3) + S_NH = cmps.S_NH4.copy('S_NH') + S_NH.description = ('Concentration of dissolved ammonium') + + # S_NO (g N/m^3) + S_NO = cmps.S_NO3.copy('S_NO') + S_NO.description = ('Concentration of dissolved nitrate/nitrite') + + # S_P (g P/m^3) + S_P = cmps.S_PO4.copy('S_P') + S_P.description = ('Concentration of dissolved phosphorus') + + # X_N_ALG (g N/m^3) + X_N_ALG = cmps.X_B_Subst.copy('X_N_ALG') + X_N_ALG.description = 'Concentration of algal cell-associated nitrogen' + X_N_ALG.measured_as = 'N' + X_N_ALG.i_C = X_N_ALG.i_P = X_N_ALG.i_COD = X_N_ALG.f_BOD5_COD = X_N_ALG.f_uBOD_COD = X_N_ALG.f_Vmass_Totmass = 0 + X_N_ALG.i_mass = 1 + + # X_P_ALG (g P/m^3) + X_P_ALG = cmps.X_B_Subst.copy('X_P_ALG') + X_P_ALG.description = 'Concentration of algal cell-associated phosphorus' + X_P_ALG.measured_as = 'P' + X_P_ALG.i_C = X_P_ALG.i_N = X_P_ALG.i_COD = X_P_ALG.f_BOD5_COD = X_P_ALG.f_uBOD_COD = X_P_ALG.f_Vmass_Totmass = 0 + X_P_ALG.i_mass = 1 + + '''added from asm2d''' + # S_N2 (g N/m^3) + S_N2 = cmps.S_N2.copy('S_N2') + S_N2.description = ('Concentration of dinitrogen') + + # S_ALK (g C/m^3) + S_ALK = cmps.S_CO3.copy('S_ALK') # measured as g C, not as mole HCO3- + S_ALK.description = ('Concentration of alkalinity') + + # S_I (g COD/m^3) + S_I = cmps.S_U_E.copy('S_I') + S_I.description = ('Concentration of inert soluble organic material') + + # X_I (g COD/m^3) + X_I = cmps.X_U_OHO_E.copy('X_I') + X_I.description = ('Concentration of inert particulate organic material') + + # X_S (g COD/m^3) + X_S = cmps.X_B_Subst.copy('X_S') + X_S.description = ('Concentration of slowly biodegradable substrates') + + # X_H (g COD/m^3) + X_H = cmps.X_OHO.copy('X_H') + X_H.description = ('Concentration of heterotrophic organisms (including denitrifer)') + + # X_AUT (g COD/m^3) + X_AUT = cmps.X_AOO.copy('X_AUT') + X_AUT.description = ('Concentration of nitrifying organisms') + + S_I.i_N = 0.01 + # S_F.i_N = 0.03 + X_I.i_N = 0.02 + X_S.i_N = 0.04 + X_H.i_N = X_AUT.i_N = 0.07 + + S_I.i_P = 0.00 + # S_F.i_P = 0.01 + X_I.i_P = 0.01 + X_S.i_P = 0.01 + X_H.i_P = X_AUT.i_P = 0.02 + + X_I.i_mass = 0.75 + X_S.i_mass = 0.75 + X_H.i_mass = X_AUT.i_mass = 0.9 + + cmps_pm2_asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, + S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, + S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, cmps.H2O]) + + cmps_pm2_asm2d.default_compile() + + if set_thermo: settings.set_thermo(cmps_pm2_asm2d) + return cmps_pm2_asm2d + +# create_pm2_asm2d_cmps() + +#%% +# ============================================================================= +# kinetic rate functions +# ============================================================================= + +# Calculation of ratio +def ratio(numerator, denominator, minimum, maximum): + return min(max(minimum, numerator / denominator), maximum) + +# Calculation of 'I_0' (for initial sensitivity analysis using calculated I_0) +def calc_irrad(t): + ''' + :param t: time [days] + :return: I_0, calculated irradiance [uE/m^2/s] + + -Assumes 14 hours of daylight + ''' + daylight_hours = 14.0 # hours + start_time = (12.0 - daylight_hours / 2) / 24.0 # days + end_time = (12.0 + daylight_hours / 2) / 24.0 # days + if t-np.floor(t) < start_time or t-np.floor(t) > end_time: + return 0 + else: + return 400.0 * (np.sin(2 * np.pi * (((t - np.floor(t)) - 5 / 24) / (14 / 24)) - np.pi / 2) + 1) / 2 + +# Calculation of 'I' from 'I_0' (Beer-Lambert) +def attenuation(light, X_TSS, a_c, b_reactor): + ''' + :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or + photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] + :param X_TSS: total biomass concentration (X_ALG + X_CH + X_LI) * i_mass [g TSS/m^3] + :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] + :parma b_reactor: thickness of reactor along light path [m] + :return: I, depth-averaged irradiance [uE/m^2/s] + ''' + if X_TSS > 0: + i_avg = (light * (1 - np.exp(-a_c * X_TSS * b_reactor))) / (a_c * X_TSS * b_reactor) + return min(i_avg, light) + else: + return light + +# Calculation of 'f_I' from 'I' (Eilers & Peeters) +def irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt): + ''' + :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param I_opt: optimal irradiance [uE/m^2/s] + :return: f_I, irradiance response function [unitless] + ''' + if X_carbon > 0: + f_I = i_avg / (i_avg + I_n * (0.25 - (5 * X_CHL/X_carbon)) * ((i_avg ** 2 / I_opt ** 2) - (2 * i_avg / I_opt) + 1)) + return min(1, max(0, f_I)) + else: + return 0 + +# Droop model +def droop(quota, subsistence_quota, exponent): + ''' + :param quota: Q_N or Q_P [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] + :return: rate [unitless] + ''' + return 1 - (subsistence_quota / quota) ** exponent + +# Monod model +def monod(substrate, half_sat_const, exponent): + ''' + :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] + :param half_sat_const: K_N or K_P [g N or g P/m^3] + :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] + :return: rate [unitless] + ''' + return (substrate / (half_sat_const + substrate)) ** exponent + +# Temperature model (Arrhenius) +def temperature(temp, arr_a, arr_e): + ''' + :param temp: temperature (will be imported from input excel file) [K] + :param arr_a: arrhenius constant (A) (Goldman et al., 1974) [unitless] + :param arr_e: arrhenius exponential constant (E/R) (Goldman et al., 1974) [K] + :return: temperature component of overall growth equation [unitless] + ''' + return arr_a * np.exp(-arr_e / temp) # Used equation from Goldman et al., 1974 + +# Photoadaptation (_p1) +def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): + ''' + :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param k_gamma: photoadaptation coefficient [unitless] + :return: photoadaptation rate [g Chl/m^3/d] + ''' + if X_carbon > 0: + return 24 * ((0.2 * i_avg / I_n) / (k_gamma + (i_avg / I_n))) *\ + (0.01 + 0.03 * ((np.log(i_avg / I_n + 0.005)) / (np.log(0.01))) - X_CHL/X_carbon) * X_carbon + else: return 0 + +# Nutrients uptake (_p2, _p3, _p4, _p5, _p6) +def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param quota: Q_N or Q_P [g N or g P/g COD] + :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] + :param uptake_rate: V_NH, V_NO or V_P [g N or g P/g COD/d] + :param half_sat_const: K_N or K_P [g N or g P/m^3] + :param maximum_quota: Q_N_max or Q_P_max [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :return: nutrient uptake rate [g N or g P/m^3/d] + ''' + return uptake_rate * monod(substrate, half_sat_const, 1) * X_ALG * \ + ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 + +# Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose growth rate (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) +def max_total_growth(X_ALG, mu_max, f_np, f_temp): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param mu_max: maximum specific growth rate [d^(-1)] + :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] + :param f_temp: temperature correction factor (between 0 and 1) [unitless] + :return: maximum total growth rate for a particular mechanism, + without considering carbon source or light inhibition (= product of shared terms in growth-related equations) [g COD/m^3/d] + ''' + return mu_max * f_np * X_ALG * f_temp + +# Split the total growth rate between three processes (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) +def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): + ''' + :param f_I: irradiance response function (calculated from 'irrad_response' method) [unitless] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] + :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] + :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] + :param K_STO: half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013) [g COD/g COD] + :return: splits the total growth rate between three processes, + growth, growth on stored carbohydrates, and growth on stored lipid (= process-specific terms) [unitless] + ''' + numerators = np.asarray([K_STO * (1 - f_I), rho * f_CH, f_LI * Y_CH / Y_LI]) + return numerators/(sum(numerators)) + +# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) +def storage_saturation(f, f_max, beta): + ''' + :param f: f_CH or f_LI [g COD/g COD] + :param f_max: f_CH_max or f_LI_max [g COD/g COD] + :param beta: beta_1 or beta_2 [unitless] + :return: part of storage equations [unitless] + ''' + return 1 - (f / f_max) ** beta + +# Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose maintenance rate (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) +def max_total_maintenance(X_ALG, m_ATP): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param m_ATP: specific maintenance rate [g ATP/g COD/d] + :return: maximum total maintenance rate for a particular mechanism + (= product of shared terms in maintenance-related equations) [g COD/m^3/d] + ''' + return m_ATP * X_ALG + +# Split the total maintenance rate between three processes (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) +def maintenance_split(f_CH, f_LI, rho, Y_CH, Y_LI, Y_X_ALG, Y_ATP, K_STO): + ''' + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] + :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] + :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] + :param Y_X_ALG: yield of carbon-accumulating phototrophic organisms, Y_X_ALG_PHO, Y_X_ALG_HET_ACE, or Y_X_ALG_HET_GLU [g COD/g COD] + :param Y_ATP: yield of ATP, Y_ATP_PHO, Y_ATP_HET_ACE, or Y_ATP_HET_GLU [g ATP/g COD] + :param K_STO: half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013) [g COD/g COD] + :return: splits the total maintenance rate between three processes, + stored carbohydrate degradation, stored lipid degradation, and endogenous respiration (= process-specific terms) [unitless] + ''' + numerators = np.asarray([rho * f_CH, f_LI * Y_CH / Y_LI, K_STO]) + yield_ratios = np.asarray([Y_CH, Y_LI, Y_X_ALG]) / Y_ATP + return numerators/(sum(numerators)) * yield_ratios + +# Storage of carbohydrate/lipid (_p8, _p9, _p16, _p17, _p24, _p25) +def storage(X_ALG, f_np, response, saturation, storage_rate): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] + :param response: f_I (irradiance response function, calculated from 'irrad_response' method), acetate_response (monod(S_A, K_A, 1)), or glucose_response (monod(S_F, K_F, 1)) [unitless] + :param saturation: 1 - (f / f_max) ** beta (calculated from 'storage_saturation' method) [unitless] + :param storage_rate: q_CH or q_LI [g COD/g COD/d] + :return: storage rate [g COD/m^3/d] + ''' + return storage_rate * saturation * (1 - f_np) * response * X_ALG + +def rhos_pm2_asm2d(state_arr, params): + + # extract values of state variables + c_arr = state_arr[:21] + temp = state_arr[22] + light = state_arr[23] # imported from input file assumed + + # Q = state_arr[14] # Flow rate + # t = state_arr[15] # time + + # light = calc_irrad(t) # when to use calculated light (I_0) + + # X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, H2O = c_arr + X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr + + # extract values of parameters + cmps = params['cmps'] + a_c = params['a_c'] + I_n = params['I_n'] + arr_a = params['arr_a'] + arr_e = params['arr_e'] + beta_1 = params['beta_1'] + beta_2 = params['beta_2'] + b_reactor = params['b_reactor'] + I_opt = params['I_opt'] + k_gamma = params['k_gamma'] + K_N = params['K_N'] + K_P = params['K_P'] + K_A = params['K_A'] + K_F = params['K_F'] + rho = params['rho'] + K_STO = params['K_STO'] + f_CH_max = params['f_CH_max'] + f_LI_max = params['f_LI_max'] + m_ATP = params['m_ATP'] + mu_max = params['mu_max'] + q_CH = params['q_CH'] + q_LI = params['q_LI'] + Q_N_max = params['Q_N_max'] + Q_N_min = params['Q_N_min'] + Q_P_max = params['Q_P_max'] + Q_P_min = params['Q_P_min'] + V_NH = params['V_NH'] + V_NO = params['V_NO'] + V_P = params['V_P'] + exponent = params['exponent'] + Y_ATP_PHO = params['Y_ATP_PHO'] + Y_CH_PHO = params['Y_CH_PHO'] + Y_LI_PHO = params['Y_LI_PHO'] + Y_X_ALG_PHO = params['Y_X_ALG_PHO'] + Y_ATP_HET_ACE = params['Y_ATP_HET_ACE'] + Y_CH_NR_HET_ACE = params['Y_CH_NR_HET_ACE'] + Y_LI_NR_HET_ACE = params['Y_LI_NR_HET_ACE'] + Y_X_ALG_HET_ACE = params['Y_X_ALG_HET_ACE'] + Y_ATP_HET_GLU = params['Y_ATP_HET_GLU'] + Y_CH_NR_HET_GLU = params['Y_CH_NR_HET_GLU'] + Y_LI_NR_HET_GLU = params['Y_LI_NR_HET_GLU'] + Y_X_ALG_HET_GLU = params['Y_X_ALG_HET_GLU'] + n_dark = params['n_dark'] + + '''added from asm2d''' + f_SI = params['f_SI'] + Y_H = params['Y_H'] + f_XI_H = params['f_XI_H'] + Y_A = params['Y_A'] + f_XI_AUT = params['f_XI_AUT'] + K_h = params['K_h'] + eta_NO3 = params['eta_NO3'] + eta_fe = params['eta_fe'] + K_O2 = params['K_O2'] + K_NO3 = params['K_NO3'] + K_X = params['K_X'] + mu_H = params['mu_H'] + q_fe = params['q_fe'] + eta_NO3_H = params['eta_NO3_H'] + b_H = params['b_H'] + K_O2_H = params['K_O2_H'] + K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H + K_fe = params['K_fe'] + K_A_H = params['K_A_H'] + K_NO3_H = params['K_NO3_H'] + K_NH4_H = params['K_NH4_H'] + K_P_H = params['K_P_H'] + K_ALK_H = params['K_ALK_H'] + mu_AUT = params['mu_AUT'] + b_AUT = params['b_AUT'] + K_O2_AUT = params['K_O2_AUT'] + K_NH4_AUT = params['K_NH4_AUT'] + K_ALK_AUT = params['K_ALK_AUT'] + K_P_AUT = params['K_P_AUT'] + +# intermediate variables + f_CH = ratio(X_CH, X_ALG, 0, f_CH_max) + f_LI = ratio(X_LI, X_ALG, 0, f_LI_max) + + # Q_N = ratio(X_N_ALG, X_ALG, Q_N_min, Q_N_max) + # Q_P = ratio(X_P_ALG, X_ALG, Q_P_min, Q_P_max) + + alg_iN, alg_iP = cmps.X_ALG.i_N, cmps.X_ALG.i_P + Q_N = ratio(X_N_ALG+X_ALG*alg_iN, X_ALG, Q_N_min, Q_N_max) + Q_P = ratio(X_P_ALG+X_ALG*alg_iP, X_ALG, Q_P_min, Q_P_max) + + idx = cmps.indices(['X_ALG', 'X_CH', 'X_LI']) + X_bio = np.array([X_ALG, X_CH, X_LI]) + X_TSS = sum(X_bio * cmps.i_mass[idx]) + X_carbon = sum(X_bio * cmps.i_C[idx]) + + i_avg = attenuation(light, X_TSS, a_c, b_reactor) + f_I = irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt) + dark_response = max(f_I, n_dark) + acetate_response = monod(S_A, K_A, 1) + glucose_response = monod(S_F, K_F, 1) + + f_np = min(droop(Q_N, Q_N_min, exponent), droop(Q_P, Q_P_min, exponent)) + f_temp = temperature(temp, arr_a, arr_e) + + f_sat_CH = storage_saturation(f_CH, f_CH_max, beta_1) + f_sat_LI = storage_saturation(f_LI, f_LI_max, beta_2) + + max_total_growth_rho = max_total_growth(X_ALG, mu_max, f_np, f_temp) + max_maintenance_rho = max_total_maintenance(X_ALG, m_ATP) + # light = calc_irrad(t) + + # calculate kinetic rate values + rhos = np.empty(30) + + rhos[0] = photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma) + + rhos[1] = nutrient_uptake(X_ALG, Q_N, S_NH, V_NH, K_N, Q_N_max, Q_N_min) + rhos[[2,3,4]] = nutrient_uptake(X_ALG, Q_N, S_NO, V_NO, K_N, Q_N_max, Q_N_min) * (K_N/(K_N + S_NH)) + + rhos[5] = nutrient_uptake(X_ALG, Q_P, S_P, V_P, K_P, Q_P_max, Q_P_min) + + rhos[[6,9,10]] = max_total_growth_rho \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) + rhos[6] *= f_I + rhos[[9,10]] *= dark_response + + rhos[[14,17,18]] = max_total_growth_rho \ + * acetate_response \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, K_STO) + + rhos[[22,25,26]] = max_total_growth_rho \ + * glucose_response \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, K_STO) + + rhos[[11,12,13]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, Y_ATP_PHO, K_STO) + + rhos[[19,20,21]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_ATP_HET_ACE, K_STO) + + rhos[[27,28,29]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU, Y_ATP_HET_GLU, K_STO) + + rhos[7] = storage(X_ALG, f_np, f_I, f_sat_CH, q_CH) + rhos[8] = storage(X_ALG, f_np, f_I, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[15] = storage(X_ALG, f_np, acetate_response, f_sat_CH, q_CH) + rhos[16] = storage(X_ALG, f_np, acetate_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[23] = storage(X_ALG, f_np, glucose_response, f_sat_CH, q_CH) + rhos[24] = storage(X_ALG, f_np, glucose_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + + return rhos + +#%% +# ============================================================================= +# PM2 class +# ============================================================================= + +class PM2ASM2d(CompiledProcesses): + ''' + Parameters + ---------- + components: class:`CompiledComponents`, optional + Components corresponding to each entry in the stoichiometry array, + defaults to thermosteam.settings.chemicals. + a_c : float, optional + PAR absorption coefficient on a TSS (total suspended solids) basis, in [m^2/g TSS]. + The default is 0.049. + I_n : float, optional + Maximum incident PAR irradiance (“irradiance at noon”), in [uE/m^2/s]. + The default is 250. + arr_a : float, optional + Arrhenius constant (A), in [unitless]. + The default is 1.8 * 10**10. + arr_e : float, optional + Arrhenius exponential constant (E/R), in [K]. + The default is 6842. + beta_1 : float, optional + Power coefficient for carbohydrate storage inhibition, in [unitless]. + The default is 2.90. + beta_2 : float, optional + Power coefficient for lipid storage inhibition, in [unitless]. + The default is 3.50. + b_reactor : float, optional + Thickness of reactor along light path, in [m]. + The default is 0.03. + I_opt : float, optional + Optimal irradiance, in [uE/m^2/s]. + The default is 300. + k_gamma : float, optional + Photoadaptation coefficient, in [unitless]. + The default is 0.00001. + K_N : float, optional + Nitrogen half-saturation constant, in [g N/m^3]. + The default is 0.1. + K_P : float, optional + Phosphorus half-saturation constant, in [g P/m^3]. + The default is 1.0. + K_A : float, optional + Organic carbon half-saturation constant (acetate) (Wagner, 2016), in [g COD/m^3]. + The default is 6.3. + K_F : float, optional + Organic carbon half-saturation constant (glucose); assumes K_A = K_F, in [g COD/m^3]. + The default is 6.3. + rho : float, optional + Carbohydrate relative preference factor (calibrated in Guest et al., 2013), in [unitless]. + The default is 1.186. + K_STO : float, optional + Half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013), in [g COD/g COD]. + The default is 1.566. + f_CH_max : float, optional + Maximum achievable ratio of stored carbohydrates to functional cells, in [g COD/g COD]. + The default is 0.819. + f_LI_max : float, optional + Maximum achievable ratio of stored lipids to functional cells, in [g COD/g COD]. + The default is 3.249. + m_ATP : float, optional + Specific maintenance rate, in [g ATP/g COD/d]. + The default is 15.835. + mu_max : float, optional + Maximum specific growth rate, in [d^(-1)]. + The default is 1.969. + q_CH : float, optional + Maximum specific carbohydrate storage rate, in [g COD/g COD/d]. + The default is 0.594. + q_LI : float, optional + Maximum specific lipid storage rate, in [g COD/g COD/d]. + The default is 0.910. + Q_N_max : float, optional + Maximum nitrogen quota, in [g N/g COD]. + The default is 0.417. + Q_N_min : float, optional + Nitrogen subsistence quota, in [g N/g COD]. + The default is 0.082. + Q_P_max : float, optional + Maximum phosphorus quota, in [g P/g COD]. + The default is 0.092. + Q_P_min : float, optional + Phosphorus subsistence quota; assumes N:P ratio of 5:1, in [g P/g COD]. + The default is 0.0163. + V_NH : float, optional + Maximum specific ammonium uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. + The default is 0.254. + V_NO : float, optional + Maximum specific nitrate uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. + The default is 0.254. + V_P : float, optional + Maximum specific phosphorus uptake rate (calibrated in Guest et al., 2013), in [g P/g COD/d]. + The default is 0.016. + exponent : float, optional + Exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013), in [unitless] + The default is 4. + Y_ATP_PHO : float, optional + Yield of ATP on CO2 fixed to G3P, in [g ATP/g CO2]. + The default is 55.073. + Y_CH_PHO : float, optional + Yield of storage carbohydrate (as polyglucose, PG) on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.754. + Y_LI_PHO : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.901. + Y_X_ALG_PHO : float, optional + Yield of carbon-accumulating phototrophic organisms on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.450. + Y_ATP_HET_ACE : float, optional + Yield of ATP on acetate fixed to acetyl-CoA, in [g ATP/g COD]. + The default is 39.623. + Y_CH_NR_HET_ACE : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on acetate fixed to acetyl-CoA under nutrient-replete condition, in [g COD/g COD]. + The default is 0.625. + Y_CH_ND_HET_ACE : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on acetate fixed to acetyl-CoA under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.600. + Y_LI_NR_HET_ACE : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on acetate fixed to acetyl-CoA under nutrient-replete condition, in [g COD/g COD]. + The default is 1.105. + Y_LI_ND_HET_ACE : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on acetate fixed to acetyl-CoA under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.713. + Y_X_ALG_HET_ACE : float, optional + Yield of carbon-accumulating phototrophic organisms on acetate fixed to acetyl-CoA, in [g COD/g COD]. + The default is 0.216. + Y_ATP_HET_GLU : float, optional + Yield of ATP on glucose fixed to G6P, in [g ATP/g COD]. + The default is 58.114. + Y_CH_NR_HET_GLU : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on glucose fixed to G6P under nutrient-replete condition, in [g COD/g COD]. + The default is 0.917. + Y_CH_ND_HET_GLU : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on glucose fixed to G6P under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.880. + Y_LI_NR_HET_GLU : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on glucose fixed to G6P under nutrient-replete condition, in [g COD/g COD]. + The default is 1.620. + Y_LI_ND_HET_GLU : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on glucose fixed to G6P under nutrient-deplete condition, in [g COD/g COD]. + The default is 1.046. + Y_X_ALG_HET_GLU : float, optional + Yield of carbon-accumulating phototrophic organisms on glucose fixed to G6P, in [g COD/g COD]. + The default is 0.317. + n_dark: float, optional + Dark growth reduction factor, in [unitless] + The default is 0.7. + path : str, optional + Alternative file path for the Petersen matrix. + The default is None. + + Examples + -------- + >>> from qsdsan import processes as pc + >>> cmps = pc.create_pm2_cmps() + >>> pm2 = pc.PM2() + >>> pm2.show() + PM2([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, + growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, + carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, + growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, + carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, + growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, + carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu]) + ''' + + _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', + 'Y_CH_NR_HET_ACE', 'Y_LI_NR_HET_ACE', 'Y_X_ALG_HET_ACE', + 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') + + _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', + *_shared_params) + + _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', + 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', + 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', + 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') + + def __new__(cls, components=None, + a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, + K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, + f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, + Q_N_max=0.417, Q_N_min=0.082, Q_P_max=0.092, Q_P_min=0.0163, V_NH=0.254, V_NO=0.254, V_P=0.016, exponent=4, + Y_ATP_PHO=55.073, Y_CH_PHO=0.754, Y_LI_PHO=0.901, Y_X_ALG_PHO=0.450, + Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, + Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, + Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, + Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, + path=None, **kwargs): + + if not path: path = _path + self = Processes.load_from_file(path, + components=components, + conserved_for=('COD', 'C', 'N', 'P'), + parameters=cls._stoichio_params, + compile=False) + + if path == _path: + _p3 = Process('nitrate_uptake_pho', + 'S_NO -> [?]S_O2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + _p4 = Process('nitrate_uptake_ace', + 'S_NO + [?]S_A -> [?]S_CO2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + _p5 = Process('nitrate_uptake_glu', + 'S_NO + [?]S_F -> [?]S_CO2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + self.insert(2, _p3) + self.insert(3, _p4) + self.insert(4, _p5) + + self.compile(to_class=cls) + + self.set_rate_function(rhos_pm2_asm2d) + shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, + Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, + Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) + stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, + *shared_values) + Q_N_min = max(self.Th_Q_N_min, Q_N_min) + Q_P_min = max(self.Th_Q_P_min, Q_P_min) + kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, + K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, + q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, + Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, + *shared_values, n_dark, self._components) + + dct = self.__dict__ + dct.update(kwargs) + dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_values)) + self.rate_function._params = dict(zip(cls._kinetic_params, kinetic_values)) + + return self + + def set_parameters(self, **parameters): + '''Set values to stoichiometric and/or kinetic parameters.''' + stoichio_only = {k:v for k,v in parameters.items() if k in self._stoichio_params} + self._parameters.update(stoichio_only) + if self._stoichio_lambdified is not None: + self.__dict__['_stoichio_lambdified'] = None + if 'Q_N_min' in parameters.keys(): + if parameters['Q_N_min'] < self.Th_Q_N_min: + raise ValueError(f'Value for Q_N_min must not be less than the ' + f'theoretical minimum {self.Th_Q_N_min}') + if 'Q_P_min' in parameters.keys(): + if parameters['Q_P_min'] < self.Th_Q_P_min: + raise ValueError(f'Value for Q_P_min must not be less than the ' + f'theoretical minimum {self.Th_Q_P_min}') + self.rate_function.set_param(**parameters) + + @property + def Th_Q_N_min(self): + return abs(self.stoichiometry.loc['growth_pho', 'X_N_ALG'])*1.001 + + @property + def Th_Q_P_min(self): + return abs(self.stoichiometry.loc['growth_pho', 'X_P_ALG'])*1.001 \ No newline at end of file From d9147984eaec8cd1df1938182398db2f406489ee Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim Date: Mon, 22 Apr 2024 22:26:31 -0500 Subject: [PATCH 02/22] added rate equations rate equations updated final class underway --- .../{_pm2_asm2d.tsv => _pm2asm2d.tsv} | 0 .../processes/{_pm2_asm2d.py => _pm2asm2d.py} | 404 +++++++++++------- 2 files changed, 243 insertions(+), 161 deletions(-) rename qsdsan/data/process_data/{_pm2_asm2d.tsv => _pm2asm2d.tsv} (100%) rename qsdsan/processes/{_pm2_asm2d.py => _pm2asm2d.py} (79%) diff --git a/qsdsan/data/process_data/_pm2_asm2d.tsv b/qsdsan/data/process_data/_pm2asm2d.tsv similarity index 100% rename from qsdsan/data/process_data/_pm2_asm2d.tsv rename to qsdsan/data/process_data/_pm2asm2d.tsv diff --git a/qsdsan/processes/_pm2_asm2d.py b/qsdsan/processes/_pm2asm2d.py similarity index 79% rename from qsdsan/processes/_pm2_asm2d.py rename to qsdsan/processes/_pm2asm2d.py index 08a0a294..7a0d98b2 100644 --- a/qsdsan/processes/_pm2_asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -16,97 +16,97 @@ from qsdsan.utils import ospath, data_path import numpy as np -__all__ = ('create_pm2_asm2d_cmps', 'PM2ASM2d') +__all__ = ('create_pm2asm2d_cmps', 'PM2ASM2d') -_path = ospath.join(data_path, 'process_data/_pm2_asm2d.tsv') +_path = ospath.join(data_path, 'process_data/_pm2asm2d.tsv') # _load_components = settings.get_default_chemicals -#%% +#%% # ============================================================================= # PM2ASM2d-specific components # ============================================================================= -def create_pm2_asm2d_cmps(set_thermo=True): +def create_pm2asm2d_cmps(set_thermo=True): cmps = Components.load_default() # X_CHL (g Chl/m^3) X_CHL = Component(ID = 'X_CHL', - formula = 'C55H72MgN4O5', + formula = 'C55H72MgN4O5', description = 'Chlorophyll content of cells', - particle_size = 'Particulate', - degradability = 'Slowly', - organic = True) - + particle_size = 'Particulate', + degradability = 'Slowly', + organic = True) + # X_ALG (g COD/m^3) - X_ALG = cmps.X_OHO.copy('X_ALG') + X_ALG = cmps.X_OHO.copy('X_ALG') X_ALG.description = 'Concentration of carbon-accumulating mixotrophic organisms' - X_ALG.formula = 'CH1.8O0.5N0.2P0.018' + X_ALG.formula = 'CH1.8O0.5N0.2P0.018' X_ALG.f_BOD5_COD = X_ALG.f_uBOD_COD = None - X_ALG.f_Vmass_Totmass = 0.89 - + X_ALG.f_Vmass_Totmass = 0.89 + # X_CH (g COD/m^3) - X_CH = cmps.X_GAO_Gly.copy('X_CH') + X_CH = cmps.X_GAO_Gly.copy('X_CH') X_CH.description = 'Concentration of stored carbohydrates' X_CH.formula = 'CH2O' X_CH.f_BOD5_COD = X_CH.f_uBOD_COD = None - + # X_LI (g COD/m^3) - X_LI = cmps.X_GAO_Gly.copy('X_LI') + X_LI = cmps.X_GAO_Gly.copy('X_LI') X_LI.description = 'Concentration of stored lipids' X_LI.formula = 'CH1.92O0.118' X_LI.f_BOD5_COD = X_LI.f_uBOD_COD = None - + # S_CO2 (g CO2/m^3) - S_CO2 = Component.from_chemical(ID = 'S_CO2', + S_CO2 = Component.from_chemical(ID = 'S_CO2', chemical = 'CO2', description = 'Soluble carbon dioxide', particle_size = 'Soluble', - degradability = 'Undegradable', - organic = False) - + degradability = 'Undegradable', + organic = False) + # S_A (g COD/m^3) - S_A = cmps.S_Ac.copy('S_A') + S_A = cmps.S_Ac.copy('S_A') S_A.description = 'Concentration of extracellular dissolved organic carbon (acetate)' - + # S_F (g COD/m^3) S_F = Component.from_chemical(ID = 'S_F', - chemical = 'glucose', + chemical = 'glucose', description = 'Concentration of extracellular dissolved organic carbon (glucose)', - measured_as = 'COD', + measured_as = 'COD', particle_size = 'Soluble', - degradability = 'Readily', - organic = True) - + degradability = 'Readily', + organic = True) + # S_O2 (g O2/m^3) - S_O2 = cmps.S_O2.copy('S_O2') - S_O2.description = ('Concentration of dissolved oxygen') - + S_O2 = cmps.S_O2.copy('S_O2') + S_O2.description = ('Concentration of dissolved oxygen') + # S_NH (g N/m^3) - S_NH = cmps.S_NH4.copy('S_NH') - S_NH.description = ('Concentration of dissolved ammonium') - + S_NH = cmps.S_NH4.copy('S_NH') + S_NH.description = ('Concentration of dissolved ammonium') + # S_NO (g N/m^3) - S_NO = cmps.S_NO3.copy('S_NO') + S_NO = cmps.S_NO3.copy('S_NO') S_NO.description = ('Concentration of dissolved nitrate/nitrite') - + # S_P (g P/m^3) - S_P = cmps.S_PO4.copy('S_P') + S_P = cmps.S_PO4.copy('S_P') S_P.description = ('Concentration of dissolved phosphorus') - + # X_N_ALG (g N/m^3) X_N_ALG = cmps.X_B_Subst.copy('X_N_ALG') X_N_ALG.description = 'Concentration of algal cell-associated nitrogen' X_N_ALG.measured_as = 'N' X_N_ALG.i_C = X_N_ALG.i_P = X_N_ALG.i_COD = X_N_ALG.f_BOD5_COD = X_N_ALG.f_uBOD_COD = X_N_ALG.f_Vmass_Totmass = 0 X_N_ALG.i_mass = 1 - + # X_P_ALG (g P/m^3) X_P_ALG = cmps.X_B_Subst.copy('X_P_ALG') X_P_ALG.description = 'Concentration of algal cell-associated phosphorus' X_P_ALG.measured_as = 'P' X_P_ALG.i_C = X_P_ALG.i_N = X_P_ALG.i_COD = X_P_ALG.f_BOD5_COD = X_P_ALG.f_uBOD_COD = X_P_ALG.f_Vmass_Totmass = 0 X_P_ALG.i_mass = 1 - + '''added from asm2d''' # S_N2 (g N/m^3) S_N2 = cmps.S_N2.copy('S_N2') @@ -151,18 +151,18 @@ def create_pm2_asm2d_cmps(set_thermo=True): X_I.i_mass = 0.75 X_S.i_mass = 0.75 X_H.i_mass = X_AUT.i_mass = 0.9 - - cmps_pm2_asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, - S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, + + cmps_pm2asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, + S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, cmps.H2O]) - - cmps_pm2_asm2d.default_compile() - if set_thermo: settings.set_thermo(cmps_pm2_asm2d) - return cmps_pm2_asm2d + cmps_pm2asm2d.default_compile() + + if set_thermo: settings.set_thermo(cmps_pm2asm2d) + return cmps_pm2asm2d + +# create_pm2asm2d_cmps() -# create_pm2_asm2d_cmps() - #%% # ============================================================================= # kinetic rate functions @@ -191,11 +191,11 @@ def calc_irrad(t): # Calculation of 'I' from 'I_0' (Beer-Lambert) def attenuation(light, X_TSS, a_c, b_reactor): ''' - :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or - photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] + :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or + photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] :param X_TSS: total biomass concentration (X_ALG + X_CH + X_LI) * i_mass [g TSS/m^3] - :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] - :parma b_reactor: thickness of reactor along light path [m] + :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] + :parma b_reactor: thickness of reactor along light path [m] :return: I, depth-averaged irradiance [uE/m^2/s] ''' if X_TSS > 0: @@ -208,9 +208,9 @@ def attenuation(light, X_TSS, a_c, b_reactor): def irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt): ''' :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] - :param X_CHL: chlorophyll content of cells [g Chl/m^3] - :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] - :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] :param I_opt: optimal irradiance [uE/m^2/s] :return: f_I, irradiance response function [unitless] ''' @@ -227,7 +227,7 @@ def droop(quota, subsistence_quota, exponent): :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] :return: rate [unitless] - ''' + ''' return 1 - (subsistence_quota / quota) ** exponent # Monod model @@ -240,25 +240,25 @@ def monod(substrate, half_sat_const, exponent): ''' return (substrate / (half_sat_const + substrate)) ** exponent -# Temperature model (Arrhenius) +# Temperature model (Arrhenius) def temperature(temp, arr_a, arr_e): ''' - :param temp: temperature (will be imported from input excel file) [K] + :param temp: temperature (will be imported from input excel file) [K] :param arr_a: arrhenius constant (A) (Goldman et al., 1974) [unitless] :param arr_e: arrhenius exponential constant (E/R) (Goldman et al., 1974) [K] :return: temperature component of overall growth equation [unitless] ''' return arr_a * np.exp(-arr_e / temp) # Used equation from Goldman et al., 1974 -# Photoadaptation (_p1) +# Photoadaptation (_p1) def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): ''' :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] :param X_CHL: chlorophyll content of cells [g Chl/m^3] :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] - :param k_gamma: photoadaptation coefficient [unitless] - :return: photoadaptation rate [g Chl/m^3/d] + :param k_gamma: photoadaptation coefficient [unitless] + :return: photoadaptation rate [g Chl/m^3/d] ''' if X_carbon > 0: return 24 * ((0.2 * i_avg / I_n) / (k_gamma + (i_avg / I_n))) *\ @@ -266,19 +266,19 @@ def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): else: return 0 # Nutrients uptake (_p2, _p3, _p4, _p5, _p6) -def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): +def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): ''' :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] - :param quota: Q_N or Q_P [g N or g P/g COD] + :param quota: Q_N or Q_P [g N or g P/g COD] :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] :param uptake_rate: V_NH, V_NO or V_P [g N or g P/g COD/d] :param half_sat_const: K_N or K_P [g N or g P/m^3] :param maximum_quota: Q_N_max or Q_P_max [g N or g P/g COD] - :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] :return: nutrient uptake rate [g N or g P/m^3/d] ''' return uptake_rate * monod(substrate, half_sat_const, 1) * X_ALG * \ - ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 + ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 # Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose growth rate (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) def max_total_growth(X_ALG, mu_max, f_np, f_temp): @@ -287,7 +287,7 @@ def max_total_growth(X_ALG, mu_max, f_np, f_temp): :param mu_max: maximum specific growth rate [d^(-1)] :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] :param f_temp: temperature correction factor (between 0 and 1) [unitless] - :return: maximum total growth rate for a particular mechanism, + :return: maximum total growth rate for a particular mechanism, without considering carbon source or light inhibition (= product of shared terms in growth-related equations) [g COD/m^3/d] ''' return mu_max * f_np * X_ALG * f_temp @@ -296,8 +296,8 @@ def max_total_growth(X_ALG, mu_max, f_np, f_temp): def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): ''' :param f_I: irradiance response function (calculated from 'irrad_response' method) [unitless] - :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] - :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] @@ -308,7 +308,7 @@ def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): numerators = np.asarray([K_STO * (1 - f_I), rho * f_CH, f_LI * Y_CH / Y_LI]) return numerators/(sum(numerators)) -# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) +# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) def storage_saturation(f, f_max, beta): ''' :param f: f_CH or f_LI [g COD/g COD] @@ -323,7 +323,7 @@ def max_total_maintenance(X_ALG, m_ATP): ''' :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] :param m_ATP: specific maintenance rate [g ATP/g COD/d] - :return: maximum total maintenance rate for a particular mechanism + :return: maximum total maintenance rate for a particular mechanism (= product of shared terms in maintenance-related equations) [g COD/m^3/d] ''' return m_ATP * X_ALG @@ -331,8 +331,8 @@ def max_total_maintenance(X_ALG, m_ATP): # Split the total maintenance rate between three processes (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) def maintenance_split(f_CH, f_LI, rho, Y_CH, Y_LI, Y_X_ALG, Y_ATP, K_STO): ''' - :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] - :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] @@ -358,8 +358,36 @@ def storage(X_ALG, f_np, response, saturation, storage_rate): ''' return storage_rate * saturation * (1 - f_np) * response * X_ALG -def rhos_pm2_asm2d(state_arr, params): - +'''added from asm2d''' + +# Hydrolysis (_p31, _p32, _p33) +def hydrolysis(X_S, X_H, K_h, K_X): + ''' + :param X_S: concentration of slowly biodegradable substrates + :param X_H: concentration of heterotrophic organisms (including denitrifer) + :param K_h: hydrolysis rate constant + :param K_X: slowly biodegradable substrate half saturation coefficient for hydrolysis + :return: shared parts of hydrolysis equations + ''' + return K_h * (X_S/X_H) / (K_X + X_S/X_H) * X_H + +# Growth in ASM2d (_p34, _p35, _p36, _p37, _p40) +def growth_asm2d(S_NH, S_P, S_ALK, mu, X, K_NH4, K_P, K_ALK): + ''' + :param S_NH: concentration of dissolved ammonium + :param S_P: concentration of dissolved phosphorus + :param S_ALK: concentration of alkalinity + :param mu: maximum specific growth rate (mu_H or mu_AUT) + :param X: concentration of biomass (X_H or X_AUT) + :param K_NH4: ammonium (nutrient) half saturation coefficient (K_NH4_H or K_NH4_AUT) + :param K_P: phosphorus (nutrient) half saturation coefficient (K_P_H or K_P_AUT) + :param K_ALK: alkalinity half saturation coefficient (K_ALK_H or K_ALK_AUT) + :return: shared parts of growth-related equations + ''' + return mu * S_NH/(K_NH4+S_NH) * S_P/(K_P+S_P) * S_ALK/(K_ALK + S_ALK) * X + +def rhos_pm2asm2d(state_arr, params): + # extract values of state variables c_arr = state_arr[:21] temp = state_arr[22] @@ -368,26 +396,23 @@ def rhos_pm2_asm2d(state_arr, params): # Q = state_arr[14] # Flow rate # t = state_arr[15] # time - # light = calc_irrad(t) # when to use calculated light (I_0) - - # X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, H2O = c_arr X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr - + # extract values of parameters - cmps = params['cmps'] + cmps = params['cmps'] a_c = params['a_c'] I_n = params['I_n'] arr_a = params['arr_a'] arr_e = params['arr_e'] beta_1 = params['beta_1'] - beta_2 = params['beta_2'] + beta_2 = params['beta_2'] b_reactor = params['b_reactor'] I_opt = params['I_opt'] k_gamma = params['k_gamma'] K_N = params['K_N'] K_P = params['K_P'] K_A = params['K_A'] - K_F = params['K_F'] + K_F = params['K_F'] rho = params['rho'] K_STO = params['K_STO'] f_CH_max = params['f_CH_max'] @@ -417,7 +442,7 @@ def rhos_pm2_asm2d(state_arr, params): Y_LI_NR_HET_GLU = params['Y_LI_NR_HET_GLU'] Y_X_ALG_HET_GLU = params['Y_X_ALG_HET_GLU'] n_dark = params['n_dark'] - + '''added from asm2d''' f_SI = params['f_SI'] Y_H = params['Y_H'] @@ -452,17 +477,17 @@ def rhos_pm2_asm2d(state_arr, params): # intermediate variables f_CH = ratio(X_CH, X_ALG, 0, f_CH_max) f_LI = ratio(X_LI, X_ALG, 0, f_LI_max) - + # Q_N = ratio(X_N_ALG, X_ALG, Q_N_min, Q_N_max) # Q_P = ratio(X_P_ALG, X_ALG, Q_P_min, Q_P_max) - + alg_iN, alg_iP = cmps.X_ALG.i_N, cmps.X_ALG.i_P Q_N = ratio(X_N_ALG+X_ALG*alg_iN, X_ALG, Q_N_min, Q_N_max) Q_P = ratio(X_P_ALG+X_ALG*alg_iP, X_ALG, Q_P_min, Q_P_max) - + idx = cmps.indices(['X_ALG', 'X_CH', 'X_LI']) X_bio = np.array([X_ALG, X_CH, X_LI]) - X_TSS = sum(X_bio * cmps.i_mass[idx]) + X_TSS = sum(X_bio * cmps.i_mass[idx]) X_carbon = sum(X_bio * cmps.i_C[idx]) i_avg = attenuation(light, X_TSS, a_c, b_reactor) @@ -480,33 +505,33 @@ def rhos_pm2_asm2d(state_arr, params): max_total_growth_rho = max_total_growth(X_ALG, mu_max, f_np, f_temp) max_maintenance_rho = max_total_maintenance(X_ALG, m_ATP) # light = calc_irrad(t) - + # calculate kinetic rate values - rhos = np.empty(30) - + rhos = np.empty(41) + rhos[0] = photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma) - + rhos[1] = nutrient_uptake(X_ALG, Q_N, S_NH, V_NH, K_N, Q_N_max, Q_N_min) rhos[[2,3,4]] = nutrient_uptake(X_ALG, Q_N, S_NO, V_NO, K_N, Q_N_max, Q_N_min) * (K_N/(K_N + S_NH)) rhos[5] = nutrient_uptake(X_ALG, Q_P, S_P, V_P, K_P, Q_P_max, Q_P_min) - + rhos[[6,9,10]] = max_total_growth_rho \ - * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) rhos[6] *= f_I - rhos[[9,10]] *= dark_response - + rhos[[9,10]] *= dark_response + rhos[[14,17,18]] = max_total_growth_rho \ * acetate_response \ * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, K_STO) - + rhos[[22,25,26]] = max_total_growth_rho \ * glucose_response \ * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, K_STO) - + rhos[[11,12,13]] = max_maintenance_rho \ * maintenance_split(f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, Y_ATP_PHO, K_STO) - + rhos[[19,20,21]] = max_maintenance_rho \ * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_ATP_HET_ACE, K_STO) @@ -520,11 +545,25 @@ def rhos_pm2_asm2d(state_arr, params): rhos[23] = storage(X_ALG, f_np, glucose_response, f_sat_CH, q_CH) rhos[24] = storage(X_ALG, f_np, glucose_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[30] = hydrolysis(X_S, X_H, K_h, K_X) * monod(S_O2, K_O2, 1) + rhos[31] = hydrolysis(X_S, X_H, K_h, K_X) * eta_NO3 * monod(K_O2, S_O2, 1) * monod(S_NO, K_NO3, 1) + rhos[32] = hydrolysis(X_S, X_H, K_h, K_X) * eta_fe * monod(K_O2, S_O2, 1) * monod(K_NO3, S_NO, 1) + + rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[34] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) + rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[36] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) + rhos[37] = q_fe * monod(K_O2_H, S_O2, 1) * monod(K_NO3_H, S_NO, 1) * monod(S_F, K_fe, 1) * monod(S_ALK, K_ALK_H, 1) * X_H + rhos[38] = b_H * X_H + + rhos[39] = growth_asm2d(S_NH, S_P, S_ALK, mu_AUT, X_AUT, K_NH4_AUT, K_P_AUT, K_ALK_AUT) * monod(S_O2, K_O2_AUT, 1) + rhos[40] = b_AUT * X_AUT + return rhos #%% # ============================================================================= -# PM2 class +# PM2ASM2d class # ============================================================================= class PM2ASM2d(CompiledProcesses): @@ -533,82 +572,82 @@ class PM2ASM2d(CompiledProcesses): ---------- components: class:`CompiledComponents`, optional Components corresponding to each entry in the stoichiometry array, - defaults to thermosteam.settings.chemicals. + defaults to thermosteam.settings.chemicals. a_c : float, optional PAR absorption coefficient on a TSS (total suspended solids) basis, in [m^2/g TSS]. - The default is 0.049. + The default is 0.049. I_n : float, optional Maximum incident PAR irradiance (“irradiance at noon”), in [uE/m^2/s]. - The default is 250. + The default is 250. arr_a : float, optional Arrhenius constant (A), in [unitless]. - The default is 1.8 * 10**10. + The default is 1.8 * 10**10. arr_e : float, optional Arrhenius exponential constant (E/R), in [K]. - The default is 6842. + The default is 6842. beta_1 : float, optional Power coefficient for carbohydrate storage inhibition, in [unitless]. - The default is 2.90. + The default is 2.90. beta_2 : float, optional Power coefficient for lipid storage inhibition, in [unitless]. - The default is 3.50. + The default is 3.50. b_reactor : float, optional Thickness of reactor along light path, in [m]. - The default is 0.03. + The default is 0.03. I_opt : float, optional Optimal irradiance, in [uE/m^2/s]. - The default is 300. + The default is 300. k_gamma : float, optional Photoadaptation coefficient, in [unitless]. - The default is 0.00001. + The default is 0.00001. K_N : float, optional Nitrogen half-saturation constant, in [g N/m^3]. - The default is 0.1. + The default is 0.1. K_P : float, optional Phosphorus half-saturation constant, in [g P/m^3]. - The default is 1.0. + The default is 1.0. K_A : float, optional Organic carbon half-saturation constant (acetate) (Wagner, 2016), in [g COD/m^3]. - The default is 6.3. + The default is 6.3. K_F : float, optional Organic carbon half-saturation constant (glucose); assumes K_A = K_F, in [g COD/m^3]. - The default is 6.3. + The default is 6.3. rho : float, optional Carbohydrate relative preference factor (calibrated in Guest et al., 2013), in [unitless]. - The default is 1.186. + The default is 1.186. K_STO : float, optional Half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013), in [g COD/g COD]. - The default is 1.566. + The default is 1.566. f_CH_max : float, optional Maximum achievable ratio of stored carbohydrates to functional cells, in [g COD/g COD]. - The default is 0.819. + The default is 0.819. f_LI_max : float, optional Maximum achievable ratio of stored lipids to functional cells, in [g COD/g COD]. - The default is 3.249. + The default is 3.249. m_ATP : float, optional Specific maintenance rate, in [g ATP/g COD/d]. - The default is 15.835. + The default is 15.835. mu_max : float, optional Maximum specific growth rate, in [d^(-1)]. - The default is 1.969. + The default is 1.969. q_CH : float, optional Maximum specific carbohydrate storage rate, in [g COD/g COD/d]. - The default is 0.594. + The default is 0.594. q_LI : float, optional Maximum specific lipid storage rate, in [g COD/g COD/d]. - The default is 0.910. + The default is 0.910. Q_N_max : float, optional Maximum nitrogen quota, in [g N/g COD]. - The default is 0.417. + The default is 0.417. Q_N_min : float, optional Nitrogen subsistence quota, in [g N/g COD]. - The default is 0.082. + The default is 0.082. Q_P_max : float, optional Maximum phosphorus quota, in [g P/g COD]. The default is 0.092. Q_P_min : float, optional Phosphorus subsistence quota; assumes N:P ratio of 5:1, in [g P/g COD]. - The default is 0.0163. + The default is 0.0163. V_NH : float, optional Maximum specific ammonium uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. The default is 0.254. @@ -626,10 +665,10 @@ class PM2ASM2d(CompiledProcesses): The default is 55.073. Y_CH_PHO : float, optional Yield of storage carbohydrate (as polyglucose, PG) on CO2 fixed to G3P, in [g COD/g CO2]. - The default is 0.754. + The default is 0.754. Y_LI_PHO : float, optional Yield of storage lipids (as triacylglycerol, TAG) on CO2 fixed to G3P, in [g COD/g CO2]. - The default is 0.901. + The default is 0.901. Y_X_ALG_PHO : float, optional Yield of carbon-accumulating phototrophic organisms on CO2 fixed to G3P, in [g COD/g CO2]. The default is 0.450. @@ -672,46 +711,89 @@ class PM2ASM2d(CompiledProcesses): n_dark: float, optional Dark growth reduction factor, in [unitless] The default is 0.7. + + + + + f_SI = params['f_SI'] + Y_H = params['Y_H'] + f_XI_H = params['f_XI_H'] + Y_A = params['Y_A'] + f_XI_AUT = params['f_XI_AUT'] + K_h = params['K_h'] + eta_NO3 = params['eta_NO3'] + eta_fe = params['eta_fe'] + K_O2 = params['K_O2'] + K_NO3 = params['K_NO3'] + K_X = params['K_X'] + mu_H = params['mu_H'] + q_fe = params['q_fe'] + eta_NO3_H = params['eta_NO3_H'] + b_H = params['b_H'] + K_O2_H = params['K_O2_H'] + K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H + K_fe = params['K_fe'] + K_A_H = params['K_A_H'] + K_NO3_H = params['K_NO3_H'] + K_NH4_H = params['K_NH4_H'] + K_P_H = params['K_P_H'] + K_ALK_H = params['K_ALK_H'] + mu_AUT = params['mu_AUT'] + b_AUT = params['b_AUT'] + K_O2_AUT = params['K_O2_AUT'] + K_NH4_AUT = params['K_NH4_AUT'] + K_ALK_AUT = params['K_ALK_AUT'] + K_P_AUT = params['K_P_AUT'] + + + + + + + path : str, optional - Alternative file path for the Petersen matrix. + Alternative file path for the Petersen matrix. The default is None. Examples -------- >>> from qsdsan import processes as pc - >>> cmps = pc.create_pm2_cmps() - >>> pm2 = pc.PM2() - >>> pm2.show() - PM2([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, - growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, - carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, - growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, - carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, - growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, - carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu]) + >>> cmps = pc.create_pm2asm2d_cmps() + >>> pm2asm2d = pc.PM2ASM2d() + >>> pm2asm2d.show() + PM2ASM2d([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, + growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, + carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, + growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, + carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, + growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, + carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu, + aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, + hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, + auto_aero_growth, auto_lysis]) ''' _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', 'Y_CH_NR_HET_ACE', 'Y_LI_NR_HET_ACE', 'Y_X_ALG_HET_ACE', 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') - + _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', *_shared_params) - - _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', - 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', - 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', + + _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', + 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', + 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') def __new__(cls, components=None, - a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, - K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, - f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, + a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, + K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, + f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, Q_N_max=0.417, Q_N_min=0.082, Q_P_max=0.092, Q_P_min=0.0163, V_NH=0.254, V_NO=0.254, V_P=0.016, exponent=4, Y_ATP_PHO=55.073, Y_CH_PHO=0.754, Y_LI_PHO=0.901, Y_X_ALG_PHO=0.450, - Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, + Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, - Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, + Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, path=None, **kwargs): @@ -721,7 +803,7 @@ def __new__(cls, components=None, conserved_for=('COD', 'C', 'N', 'P'), parameters=cls._stoichio_params, compile=False) - + if path == _path: _p3 = Process('nitrate_uptake_pho', 'S_NO -> [?]S_O2 + X_N_ALG', @@ -740,34 +822,34 @@ def __new__(cls, components=None, components=components, ref_component='X_N_ALG', conserved_for=('COD', 'C')) - + self.insert(2, _p3) self.insert(3, _p4) self.insert(4, _p5) self.compile(to_class=cls) - - self.set_rate_function(rhos_pm2_asm2d) - shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, - Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, + + self.set_rate_function(rhos_pm2asm2d) + shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, + Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, *shared_values) Q_N_min = max(self.Th_Q_N_min, Q_N_min) Q_P_min = max(self.Th_Q_P_min, Q_P_min) - kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, - K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, - q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, + kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, + K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, + q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, *shared_values, n_dark, self._components) - + dct = self.__dict__ dct.update(kwargs) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_values)) self.rate_function._params = dict(zip(cls._kinetic_params, kinetic_values)) return self - + def set_parameters(self, **parameters): '''Set values to stoichiometric and/or kinetic parameters.''' stoichio_only = {k:v for k,v in parameters.items() if k in self._stoichio_params} @@ -775,11 +857,11 @@ def set_parameters(self, **parameters): if self._stoichio_lambdified is not None: self.__dict__['_stoichio_lambdified'] = None if 'Q_N_min' in parameters.keys(): - if parameters['Q_N_min'] < self.Th_Q_N_min: + if parameters['Q_N_min'] < self.Th_Q_N_min: raise ValueError(f'Value for Q_N_min must not be less than the ' f'theoretical minimum {self.Th_Q_N_min}') if 'Q_P_min' in parameters.keys(): - if parameters['Q_P_min'] < self.Th_Q_P_min: + if parameters['Q_P_min'] < self.Th_Q_P_min: raise ValueError(f'Value for Q_P_min must not be less than the ' f'theoretical minimum {self.Th_Q_P_min}') self.rate_function.set_param(**parameters) @@ -787,7 +869,7 @@ def set_parameters(self, **parameters): @property def Th_Q_N_min(self): return abs(self.stoichiometry.loc['growth_pho', 'X_N_ALG'])*1.001 - + @property def Th_Q_P_min(self): return abs(self.stoichiometry.loc['growth_pho', 'X_P_ALG'])*1.001 \ No newline at end of file From 293e0a6d633124823ed2beb7e2a465c9b89dfd01 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Tue, 23 Apr 2024 15:25:23 -0500 Subject: [PATCH 03/22] Update _pm2asm2d.py --- qsdsan/processes/_pm2asm2d.py | 168 ++++++++++++++++++++++++---------- 1 file changed, 121 insertions(+), 47 deletions(-) diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index 7a0d98b2..7b5adee6 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -393,8 +393,8 @@ def rhos_pm2asm2d(state_arr, params): temp = state_arr[22] light = state_arr[23] # imported from input file assumed - # Q = state_arr[14] # Flow rate - # t = state_arr[15] # time + # Q = state_arr[21] # Flow rate + # t = state_arr[22] # time X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr @@ -710,47 +710,94 @@ class PM2ASM2d(CompiledProcesses): The default is 0.317. n_dark: float, optional Dark growth reduction factor, in [unitless] - The default is 0.7. - - - - - f_SI = params['f_SI'] - Y_H = params['Y_H'] - f_XI_H = params['f_XI_H'] - Y_A = params['Y_A'] - f_XI_AUT = params['f_XI_AUT'] - K_h = params['K_h'] - eta_NO3 = params['eta_NO3'] - eta_fe = params['eta_fe'] - K_O2 = params['K_O2'] - K_NO3 = params['K_NO3'] - K_X = params['K_X'] - mu_H = params['mu_H'] - q_fe = params['q_fe'] - eta_NO3_H = params['eta_NO3_H'] - b_H = params['b_H'] - K_O2_H = params['K_O2_H'] - K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H - K_fe = params['K_fe'] - K_A_H = params['K_A_H'] - K_NO3_H = params['K_NO3_H'] - K_NH4_H = params['K_NH4_H'] - K_P_H = params['K_P_H'] - K_ALK_H = params['K_ALK_H'] - mu_AUT = params['mu_AUT'] - b_AUT = params['b_AUT'] - K_O2_AUT = params['K_O2_AUT'] - K_NH4_AUT = params['K_NH4_AUT'] - K_ALK_AUT = params['K_ALK_AUT'] - K_P_AUT = params['K_P_AUT'] - - - - - - - + The default is 0.7. + f_SI : float, optional + Production of soluble inerts in hydrolysis, in [g COD/g COD]. + The default is 0.0. + Y_H : float, optional + Heterotrophic yield coefficient, in[g COD/g COD]. + The default is 0.625. + f_XI_H : float, optional + Fraction of inert COD generated in heterotrophic biomass lysis, in [g COD/g COD]. + The default is 0.1. + Y_A : float, optional + Autotrophic yield, in [g COD/g N]. + The default is 0.24. + f_XI_AUT : float, optional + Fraction of inert COD generated in autotrophic biomass lysis, in [g COD/g COD]. + The default is 0.1. + K_h : float, optional + Hydrolysis rate constant, in [d^(-1)]. + The default is 3.0. + eta_NO3 : float, optional + Reduction factor for anoxic hydrolysis, dimensionless. + The default is 0.6. + eta_fe : float, optional + Anaerobic hydrolysis reduction factor, dimensionless. + The default is 0.4. + K_O2 : float, optional + O2 half saturation coefficient for hydrolysis, in [g O2/m^3]. + The default is 0.2. + K_NO3 : float, optional + Nitrate half saturation coefficient for hydrolysis, in [g N/m^3]. + The default is 0.5. + K_X : float, optional + Slowly biodegradable substrate half saturation coefficient for hydrolysis, in [g COD/g COD]. + The default is 0.1. + mu_H : float, optional + Heterotrophic maximum specific growth rate, in [d^(-1)]. + The default is 6.0. + q_fe : float, optional + Fermentation maximum rate, in [d^(-1)]. + The default is 3.0. + eta_NO3_H : float, optional + Reduction factor for anoxic heterotrophic growth, dimensionless. + The default is 0.8. + b_H : float, optional + Lysis and decay rate constant, in [d^(-1)]. + The default is 0.4. + K_O2_H : float, optional + O2 half saturation coefficient for heterotrophs, in [g O2/m^3]. + The default is 0.2. + K_F_H : float, optional + Fermentable substrate half saturation coefficient for heterotrophic growth (K_F in ASM2d), in [g COD/m^3]. + The default is 4.0. + K_fe : float, optional + Fermentable substrate half saturation coefficient for fermentation, in [g COD/m^3]. + The default is 4.0. + K_A_H : float, optional + VFA half saturation coefficient for heterotrophs, in [g COD/m^3]. + The default is 4.0. + K_NO3_H : float, optional + Nitrate half saturation coefficient for heterotrophs, in [g N/m^3]. + The default is 0.5. + K_NH4_H : float, optional + Ammonium (nutrient) half saturation coefficient for heterotrophs, in [g N/m^3]. + The default is 0.05. + K_P_H : float, optional + Phosphorus (nutrient) half saturation coefficient for heterotrophs, in [g P/m^3]. + The default is 0.01. + K_ALK_H : float, optional + Alkalinity half saturation coefficient for heterotrophs, in [mol(HCO3-)/m^3]. (user input unit, converted as C) + The default is 0.1. + mu_AUT : float, optional + Autotrophic maximum specific growth rate, in [d^(-1)]. + The default is 1.0. + b_AUT : float, optional + Autotrophic decay rate, in [d^(-1)]. + The default is 0.15. + K_O2_AUT : float, optional + O2 half saturation coefficient for autotrophs, in [g O2/m^3]. + The default is 0.5. + K_NH4_AUT : float, optional + Ammonium (nutrient) half saturation coefficient for autotrophs, in [g N/m^3]. + The default is 1.0. + K_ALK_AUT : float, optional + Alkalinity half saturation coefficient for autotrophs, in [mol(HCO3-)/m^3]. (user input unit, converted as C) + The default is 0.5. + K_P_AUT : float, optional + Phosphorus (nutrient) half saturation coefficient for autotrophs, in [g P/m^3]. + The default is 0.01. path : str, optional Alternative file path for the Petersen matrix. The default is None. @@ -771,6 +818,17 @@ class PM2ASM2d(CompiledProcesses): aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, auto_aero_growth, auto_lysis]) + + References + ---------- + .. [1] Henze, M.; Gujer, W.; Mino, T.; Loosdrecht, M. van. Activated Sludge + Models: ASM1, ASM2, ASM2d and ASM3; IWA task group on mathematical modelling + for design and operation of biological wastewater treatment, Ed.; IWA + Publishing: London, 2000. + .. [2] Rieger, L.; Gillot, S.; Langergraber, G.; Ohtsuki, T.; Shaw, A.; Takács, + I.; Winkler, S. Guidelines for Using Activated Sludge Models; IWA Publishing: + London, New York, 2012; Vol. 11. + https://doi.org/10.2166/9781780401164. ''' _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', @@ -778,12 +836,16 @@ class PM2ASM2d(CompiledProcesses): 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', + 'f_SI', 'Y_H', 'f_XI_H', 'Y_A', 'f_XI_AUT', *_shared_params) _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', - 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') + 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps', + 'K_h', 'eta_NO3', 'eta_fe', 'K_O2', 'K_NO3', 'K_X', 'mu_H', 'q_fe', 'eta_NO3_H', + 'b_H', 'K_O2_H', 'K_F_H', 'K_fe', 'K_A_H', 'K_NO3_H', 'K_NH4_H', 'K_P_H', + 'K_ALK_H', 'mu_AUT', 'b_AUT', 'K_O2_AUT', 'K_NH4_AUT', 'K_ALK_AUT', 'K_P_AUT') def __new__(cls, components=None, a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, @@ -795,12 +857,19 @@ def __new__(cls, components=None, Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, + f_SI=0.0, Y_H=0.625, f_XI_H=0.1, Y_A=0.24, f_XI_AUT=0.1, + K_h=3.0, eta_NO3=0.6, eta_fe=0.4, K_O2=0.2, K_NO3=0.5, K_X=0.1, + mu_H=6.0, q_fe=3.0, eta_NO3_H=0.8, b_H=0.4, K_O2_H=0.2, K_F_H=4.0, + K_fe=4.0, K_A_H=4.0, K_NO3_H=0.5, K_NH4_H=0.05, K_P_H=0.01, K_ALK_H=0.1, + mu_AUT=1.0, b_AUT=0.15, K_O2_AUT=0.5, K_NH4_AUT=1.0, K_ALK_AUT=0.5, K_P_AUT=0.01, path=None, **kwargs): if not path: path = _path + self = Processes.load_from_file(path, components=components, - conserved_for=('COD', 'C', 'N', 'P'), + conserved_for=('COD', 'C', 'N', 'P', 'charge'), + #conserved_for=('COD', 'C', 'N', 'P'), parameters=cls._stoichio_params, compile=False) @@ -834,14 +903,19 @@ def __new__(cls, components=None, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, - *shared_values) + f_SI, Y_H, f_XI_H, Y_A, f_XI_AUT, + *shared_values) Q_N_min = max(self.Th_Q_N_min, Q_N_min) Q_P_min = max(self.Th_Q_P_min, Q_P_min) kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, - *shared_values, n_dark, self._components) + *shared_values, n_dark, self._components, + K_h, eta_NO3, eta_fe, K_O2, K_NO3, K_X, mu_H, q_fe, eta_NO3_H, + b_H, K_O2_H, K_F_H, K_fe, K_A_H, K_NO3_H, K_NH4_H, K_P_H, + K_ALK_H*12, mu_AUT, b_AUT, K_O2_AUT, K_NH4_AUT, K_ALK_AUT*12, K_P_AUT, + ) dct = self.__dict__ dct.update(kwargs) From d39e31c21d45508e0a48a6a065e8aa786bd9e075 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Tue, 23 Apr 2024 17:02:10 -0500 Subject: [PATCH 04/22] split pm2 & asm2d stoichiometry path --- qsdsan/data/process_data/_pm2asm2d_1.tsv | 28 ++++++++++++++++++++++++ qsdsan/data/process_data/_pm2asm2d_2.tsv | 12 ++++++++++ qsdsan/processes/__init__.py | 4 ++++ qsdsan/processes/_pm2asm2d.py | 14 +++++++++--- 4 files changed, 55 insertions(+), 3 deletions(-) create mode 100644 qsdsan/data/process_data/_pm2asm2d_1.tsv create mode 100644 qsdsan/data/process_data/_pm2asm2d_2.tsv diff --git a/qsdsan/data/process_data/_pm2asm2d_1.tsv b/qsdsan/data/process_data/_pm2asm2d_1.tsv new file mode 100644 index 00000000..f3a977f5 --- /dev/null +++ b/qsdsan/data/process_data/_pm2asm2d_1.tsv @@ -0,0 +1,28 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +photoadaptation 1 +ammonium_uptake -1 1 +phosphorus_uptake -1 1 +growth_pho 1 ? 1 ? ? +carbohydrate_storage_pho 1 ? 1 +lipid_storage_pho 1 ? 1 +carbohydrate_growth_pho 1 (-Y_CH_PHO/Y_X_ALG_PHO) ? ? ? ? +lipid_growth_pho 1 (-Y_LI_PHO/Y_X_ALG_PHO) ? ? ? ? +carbohydrate_maintenance_pho -1 ? -1 +lipid_maintenance_pho -1 ? -1 +endogenous_respiration_pho -1 ? -1 ? ? +growth_ace 1 ? (-1)/Y_X_ALG_HET_ACE ? ? ? +carbohydrate_storage_ace 1 ? (-1)/Y_CH_ND_HET_ACE ? +lipid_storage_ace 1 ? (-1)/Y_LI_ND_HET_ACE ? +carbohydrate_growth_ace 1 (-Y_CH_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +lipid_growth_ace 1 (-Y_LI_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +carbohydrate_maintenance_ace -1 ? -1 +lipid_maintenance_ace -1 ? -1 +endogenous_respiration_ace -1 ? -1 ? ? +growth_glu 1 ? (-1)/Y_X_ALG_HET_GLU ? ? ? +carbohydrate_storage_glu 1 ? (-1)/Y_CH_ND_HET_GLU ? +lipid_storage_glu 1 ? (-1)/Y_LI_ND_HET_GLU ? +carbohydrate_growth_glu 1 (-Y_CH_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +lipid_growth_glu 1 (-Y_LI_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +carbohydrate_maintenance_glu -1 ? -1 +lipid_maintenance_glu -1 ? -1 +endogenous_respiration_glu -1 ? -1 ? ? diff --git a/qsdsan/data/process_data/_pm2asm2d_2.tsv b/qsdsan/data/process_data/_pm2asm2d_2.tsv new file mode 100644 index 00000000..8941b263 --- /dev/null +++ b/qsdsan/data/process_data/_pm2asm2d_2.tsv @@ -0,0 +1,12 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +aero_hydrolysis 1-f_SI ? ? ? f_SI -1 +anox_hydrolysis 1-f_SI ? ? ? f_SI -1 +anae_hydrolysis 1-f_SI ? ? ? f_SI -1 +hetero_growth_S_F (-1)/Y_H 1-1/Y_H ? ? ? 1 +hetero_growth_S_A (-1)/Y_H 1-1/Y_H ? ? ? 1 +denitri_S_F (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +denitri_S_A (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +ferment 1 -1 ? ? ? +hetero_lysis ? ? ? f_XI_H 1-f_XI_H -1 +auto_aero_growth (Y_A-32/7)/Y_A ? 1/Y_A ? ? 1 +auto_lysis ? ? ? f_XI_AUT 1-f_XI_AUT -1 diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index 9e0c5e9d..ddcf5d2c 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -21,6 +21,7 @@ from ._decay import * from ._kinetic_reaction import * from ._pm2 import * +from ._pm2asm2d import * from . import ( _aeration, @@ -30,6 +31,8 @@ _madm1, _decay, _kinetic_reaction, + _pm2, + _pm2asm2d, ) __all__ = ( @@ -41,4 +44,5 @@ *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, + *_pm2asm2d.__all__, ) \ No newline at end of file diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index 7b5adee6..d56cdbd8 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -18,7 +18,9 @@ __all__ = ('create_pm2asm2d_cmps', 'PM2ASM2d') -_path = ospath.join(data_path, 'process_data/_pm2asm2d.tsv') +_path = ospath.join(data_path, 'process_data/_pm2asm2d_1.tsv') +_path_2 = ospath.join(data_path, 'process_data/_pm2asm2d_2.tsv') + # _load_components = settings.get_default_chemicals #%% @@ -868,10 +870,16 @@ def __new__(cls, components=None, self = Processes.load_from_file(path, components=components, - conserved_for=('COD', 'C', 'N', 'P', 'charge'), - #conserved_for=('COD', 'C', 'N', 'P'), + conserved_for=('COD', 'C', 'N', 'P'), + parameters=cls._stoichio_params, + compile=False) + + asm2d_processes = Processes.load_from_file(_path_2, + components=components, + conserved_for=('COD', 'N', 'P', 'charge'), parameters=cls._stoichio_params, compile=False) + self.extend(asm2d_processes) if path == _path: _p3 = Process('nitrate_uptake_pho', From 4ec854adadc3a6f0ed852f03137448078ff3a0ac Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Sat, 27 Apr 2024 18:57:51 -0500 Subject: [PATCH 05/22] Update _pm2asm2d.py --- qsdsan/processes/_pm2asm2d.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index d56cdbd8..223c10bd 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -446,11 +446,11 @@ def rhos_pm2asm2d(state_arr, params): n_dark = params['n_dark'] '''added from asm2d''' - f_SI = params['f_SI'] - Y_H = params['Y_H'] - f_XI_H = params['f_XI_H'] - Y_A = params['Y_A'] - f_XI_AUT = params['f_XI_AUT'] + # f_SI = params['f_SI'] + # Y_H = params['Y_H'] + # f_XI_H = params['f_XI_H'] + # Y_A = params['Y_A'] + # f_XI_AUT = params['f_XI_AUT'] K_h = params['K_h'] eta_NO3 = params['eta_NO3'] eta_fe = params['eta_fe'] @@ -551,9 +551,9 @@ def rhos_pm2asm2d(state_arr, params): rhos[31] = hydrolysis(X_S, X_H, K_h, K_X) * eta_NO3 * monod(K_O2, S_O2, 1) * monod(S_NO, K_NO3, 1) rhos[32] = hydrolysis(X_S, X_H, K_h, K_X) * eta_fe * monod(K_O2, S_O2, 1) * monod(K_NO3, S_NO, 1) - rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F_H, 1) * monod(S_F, S_A, 1) rhos[34] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) - rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F_H, 1) * monod(S_F, S_A, 1) rhos[36] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) rhos[37] = q_fe * monod(K_O2_H, S_O2, 1) * monod(K_NO3_H, S_NO, 1) * monod(S_F, K_fe, 1) * monod(S_ALK, K_ALK_H, 1) * X_H rhos[38] = b_H * X_H From 1e1f966bec1fee221bb046d09175565583e6412f Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 12:04:17 -0700 Subject: [PATCH 06/22] 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 07/22] `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 08/22] 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 09/22] 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 10/22] 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 11/22] 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 12/22] 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 13/22] 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 14/22] 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 15/22] 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 16/22] 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 17/22] 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 18/22] 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 19/22] 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): From ccb7314a61f7c5813fc9b59a594d444edf87480e Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 14:27:33 -0500 Subject: [PATCH 20/22] temporary fix for indexer phase (thermosteam recent updates) --- qsdsan/_waste_stream.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/qsdsan/_waste_stream.py b/qsdsan/_waste_stream.py index 37c5843d..4700c430 100644 --- a/qsdsan/_waste_stream.py +++ b/qsdsan/_waste_stream.py @@ -154,7 +154,11 @@ def __init__(self, dct, F_vol, MW, phase, phase_container): def output(self, index, value): '''Concentration flows, in mg/L (g/m3).''' f_mass = value * self.MW[index] - phase = self.phase or self.phase_container.phase + if self.phase: + phase = self.phase + else: + try: phase = self.phase_container._phase + except: phase = self.phase_container if phase != 'l': raise AttributeError('Concentration only valid for liquid phase.') V_sum = self.F_vol @@ -186,7 +190,7 @@ def by_conc(self, TP): check_data=False, ) return conc -indexer.ChemicalMolarFlowIndexer.by_conc = by_conc +ChemicalMolarFlowIndexer.by_conc = by_conc del by_conc From 3f40a30e999475e0e3a1b123ebd79563d37c5000 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 14:31:10 -0500 Subject: [PATCH 21/22] add new attribute with biosteam update --- qsdsan/_sanunit.py | 1 + 1 file changed, 1 insertion(+) diff --git a/qsdsan/_sanunit.py b/qsdsan/_sanunit.py index 51cac1a7..4788b965 100644 --- a/qsdsan/_sanunit.py +++ b/qsdsan/_sanunit.py @@ -217,6 +217,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream self._assert_compatible_property_package() self._utility_cost = None + self._recycle_system = None ##### qsdsan-specific ##### for i in (*construction, *transportation, *equipment): From d2edfe64656a29e3a83acd374588f7b0d71a9ad9 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 15:39:13 -0500 Subject: [PATCH 22/22] update version --- docs/source/conf.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index cc735083..20cc44a7 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -31,7 +31,7 @@ # built documents. # # The short X.Y version. -version = '1.4.0' +version = '1.4.1' # The full version, including alpha/beta/rc tags. release = version diff --git a/setup.py b/setup.py index 2059114b..995ef6d5 100644 --- a/setup.py +++ b/setup.py @@ -18,7 +18,7 @@ setup( name='qsdsan', packages=['qsdsan'], - version='1.4.0', + version='1.4.1', license='UIUC', author='Quantitative Sustainable Design Group', author_email='quantitative.sustainable.design@gmail.com',