From 6a549ea950a004e420e6f7e2a8cf47914cc7faa8 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 12 Oct 2023 09:43:38 -0700 Subject: [PATCH] speed up temperature correction in ADM1 --- qsdsan/processes/_adm1.py | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/qsdsan/processes/_adm1.py b/qsdsan/processes/_adm1.py index cb876244..161f9cfb 100644 --- a/qsdsan/processes/_adm1.py +++ b/qsdsan/processes/_adm1.py @@ -190,6 +190,8 @@ def mass2mol_conversion(cmps): # return np.exp(theta * (T2-T1)) def T_correction_factor(T1, T2, delta_H): + """compute temperature correction factor for equilibrium constants based on + the Van't Holf equation.""" if T1 == T2: return 1 return np.exp(delta_H/(R*100) * (1/T1 - 1/T2)) # R converted to SI @@ -268,10 +270,24 @@ def rhos_adm1(state_arr, params): weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]] T_op = state_arr[-1] + if T_op == T_base: + Ka = Kab + KH = KHb / unit_conversion[7:10] + else: + T_temp = params.pop('T_op', None) + if T_op == T_temp: + params['T_op'] = T_op + Ka = params['Ka'] + KH = params['KH'] + else: + params['T_op'] = T_op + Ka = params['Ka'] = Kab * T_correction_factor(T_base, T_op, Ka_dH) + KH = params['KH'] = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] + biogas_S = state_arr[7:10].copy() biogas_p = R * T_op * state_arr[27:30] - Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) - KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] + # Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) + # KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] rhos[:-3] = ks * Cs Monod = substr_inhibit(substrates, Ks) @@ -280,12 +296,12 @@ def rhos_adm1(state_arr, params): if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu) h = brenth(acid_base_rxn, 1e-14, 1.0, - args=(weak_acids, Kas), + args=(weak_acids, Ka), xtol=1e-12, maxiter=100) # h = 10**(-7.46) - nh3 = Kas[1] * weak_acids[2] / (Kas[1] + h) - co2 = weak_acids[3] - Kas[2] * weak_acids[3] / (Kas[2] + h) + nh3 = Ka[1] * weak_acids[2] / (Ka[1] + h) + co2 = weak_acids[3] - Ka[2] * weak_acids[3] / (Ka[2] + h) biogas_S[-1] = co2 / unit_conversion[9] Iph = Hill_inhibit(h, pH_ULs, pH_LLs)