consistent MMP model between ASM and ADM

QSD-Group · Jun 25, 2024 · d41a08c · d41a08c
1 parent abc0145
commit d41a08c
Show file tree

Hide file tree

Showing 2 changed files with 36 additions and 16 deletions.
diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py
@@ -737,8 +737,8 @@ def _rhos_adm1p(state_arr, params, h=None):
     k_mmp = params['k_mmp']
     Ksp = params['Ksp']
     # K_dis = params['K_dis']
-    # K_AlOH = params['K_AlOH']
-    # K_FeOH = params['K_FeOH']
+    K_AlOH = params['K_AlOH']
+    K_FeOH = params['K_FeOH']
     S_Mg, S_Ca, X_CaCO3, X_struv, X_newb, X_ACP, X_MgCO3 = state_arr[28:35]
     X_AlOH, X_FeOH = state_arr[[35,37]]
     # f_dis = Monod(state_arr[30:35], K_dis[:5])
@@ -753,8 +753,8 @@ def _rhos_adm1p(state_arr, params, h=None):
     # if X_MgCO3 > 0: rhos_p[29] = (S_Mg * co3 - Ksp[4]) * f_dis[4]
     # else: rhos_p[29] = S_Mg * co3
 
-    # rhos_p[30] = X_AlOH * po4 * Monod(X_AlOH, K_AlOH)
-    # rhos_p[31] = X_FeOH * po4 * Monod(X_FeOH, K_FeOH)
+    rhos_p[30] = X_AlOH * po4 * Monod(X_AlOH, K_AlOH)
+    rhos_p[31] = X_FeOH * po4 * Monod(X_FeOH, K_FeOH)
 
     SI = (S_Ca * co3 / Ksp[0])**(1/2)
     if SI > 1: rhos_p[25] = X_CaCO3 * (SI-1)**2

diff --git a/qsdsan/processes/_asm2d.py b/qsdsan/processes/_asm2d.py
@@ -584,20 +584,40 @@ def _rhos_masm2d(state_arr, params, acceptor_dependent_decay=True):
     ########## precipitation-dissolution #############
     k_mmp = params['k_mmp']
     Ksp = params['Ksp']
-    K_dis = params['K_dis']
+    # K_dis = params['K_dis']
     K_AlOH = params['K_AlOH']
     K_FeOH = params['K_FeOH']
-    f_dis = Monod(state_arr[19:24], K_dis[:5])
-    if X_CaCO3 > 0: rhos[19] = (S_Ca * co3 - Ksp[0]) * f_dis[0]
-    else: rhos[19] = S_Ca * co3
-    if X_struv > 0: rhos[20] = (S_Mg * nh4 * po4 - Ksp[1]) * f_dis[1]
-    else: rhos[20] = S_Mg * nh4 * po4
-    if X_newb > 0: rhos[21] = (S_Mg * hpo4 - Ksp[2]) * f_dis[2]
-    else: rhos[21] = S_Mg * hpo4
-    if X_ACP > 0: rhos[22] = (S_Ca**3 * po4**2 - Ksp[3]) * f_dis[3]
-    else: rhos[22] = S_Ca**3 * po4**2
-    if X_MgCO3 > 0: rhos[23] = (S_Mg * co3 - Ksp[4]) * f_dis[4]
-    else: rhos[23] = S_Mg * co3
+    # f_dis = Monod(state_arr[19:24], K_dis[:5])
+    # if X_CaCO3 > 0: rhos[19] = (S_Ca * co3 - Ksp[0]) * f_dis[0]
+    # else: rhos[19] = S_Ca * co3
+    # if X_struv > 0: rhos[20] = (S_Mg * nh4 * po4 - Ksp[1]) * f_dis[1]
+    # else: rhos[20] = S_Mg * nh4 * po4
+    # if X_newb > 0: rhos[21] = (S_Mg * hpo4 - Ksp[2]) * f_dis[2]
+    # else: rhos[21] = S_Mg * hpo4
+    # if X_ACP > 0: rhos[22] = (S_Ca**3 * po4**2 - Ksp[3]) * f_dis[3]
+    # else: rhos[22] = S_Ca**3 * po4**2
+    # if X_MgCO3 > 0: rhos[23] = (S_Mg * co3 - Ksp[4]) * f_dis[4]
+    # else: rhos[23] = S_Mg * co3
+    SI = (S_Ca * co3 / Ksp[0])**(1/2)
+    if SI > 1: rhos[19] = X_CaCO3 * (SI-1)**2
+    else: rhos[19] = 0
+
+    SI = (S_Mg * nh4 * po4 / Ksp[1])**(1/3)
+    if SI > 1: rhos[20] = X_struv * (SI-1)**3
+    else: rhos[20] = 0
+
+    SI = (S_Mg * hpo4 / Ksp[2])**(1/2)
+    if SI > 1: rhos[21] =  X_newb * (SI-1)**2
+    else: rhos[21] = 0
+
+    SI = (S_Ca**3 * po4**2 / Ksp[3])**(1/5)
+    if SI > 1: rhos[22] = X_ACP * (SI-1)**5
+    else: rhos[22] = 0
+
+    SI = (S_Mg * co3 / Ksp[4])**(1/2)
+    if SI > 1: rhos[23] = X_MgCO3 * (SI-1)**2
+    else: rhos[23] = 0
+
     rhos[24] = X_AlOH * po4 * Monod(X_AlOH, K_AlOH)
     rhos[25] = X_FeOH * po4 * Monod(X_FeOH, K_FeOH)
     rhos[19:26] *= k_mmp