From ef77253e4f0c27037cf5f2ee241a2d3dac9ff8a7 Mon Sep 17 00:00:00 2001
From: Carolina Tristan <ctristan@purdue.edu>
Date: Wed, 18 Dec 2024 07:46:47 -0500
Subject: [PATCH] Add dynamic viscosity parameter to `REDstack` model and
 update calculations

---
 gdplib/reverse_electrodialysis/REDstack.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/gdplib/reverse_electrodialysis/REDstack.py b/gdplib/reverse_electrodialysis/REDstack.py
index d2055ba..7721cae 100644
--- a/gdplib/reverse_electrodialysis/REDstack.py
+++ b/gdplib/reverse_electrodialysis/REDstack.py
@@ -86,6 +86,10 @@ def build_REDstack():
     )  # Mehdizadeh, et al. (2019) Membranes, 9(6), 73. https://doi.org/10.3390/membranes9060073
     # Linear temperature dependence of the solution conductivity. The temperature coefficient of the solution conductivity is 0.02 K-1.
 
+    m.dynamic_viscosity = pyo.Param(
+        doc='Dynamic viscosity of the solution [Pa s]', default=1e-3, initialize=1e-3
+    )
+
     m.pump_eff = pyo.Param(doc='Pump efficiency [-]', default=0.75, initialize=0.75)
 
     # =============================================================================
@@ -518,7 +522,7 @@ def _pressure_x(m, x, sol):
         """
         delta_p = pyo.value(
             48
-            * ureg.convert(1, 'cP', 'Pa*s')
+            * m.dynamic_viscosity
             * ureg.convert(m.vel_init[sol], 'cm', 'm')
             / m.dh[sol] ** 2
             * m.L
@@ -548,7 +552,7 @@ def _pressure_x_b(m, x, sol):
         """
         delta_p = pyo.value(
             48
-            * ureg.convert(1, 'cP', 'Pa*s')
+            * m.dynamic_viscosity
             * ureg.convert(m.vel_ub[sol], 'cm', 'm')
             / m.dh[sol] ** 2
             * m.L
@@ -587,7 +591,7 @@ def _pressure_x_b(m, x, sol):
         bounds=lambda _, x, sol: (
             ureg.convert(
                 -48
-                * ureg.convert(1, 'cP', 'Pa*s')
+                * m.dynamic_viscosity
                 * ureg.convert(m.vel_ub[sol], 'cm', 'm')
                 / m.dh[sol] ** 2
                 * m.L,
@@ -601,7 +605,7 @@ def _pressure_x_b(m, x, sol):
             if x == m.length_domain.first()
             else ureg.convert(
                 -48
-                * ureg.convert(1, 'cP', 'Pa*s')
+                * m.dynamic_viscosity
                 * ureg.convert(m.vel_init[sol], 'cm', 'm')
                 / m.dh[sol] ** 2
                 * m.L,
@@ -977,7 +981,9 @@ def _Rcpx(m, x):
     m.PP = pyo.Var(
         domain=pyo.NonNegativeReals,
         initialize=sum(
-            ureg.convert(48e-7 * m.vel_init[sol] / m.dh[sol] ** 2, 'mbar', 'Pa')
+            48
+            * m.dynamic_viscosity
+            * ureg.convert(48e-7 * m.vel_init[sol] / m.dh[sol] ** 2, 'mbar', 'Pa')
             * m.flow_vol['rm', sol]
             / 3.6e3
             / m.pump_eff
@@ -987,7 +993,7 @@ def _Rcpx(m, x):
             None,
             sum(
                 48
-                * ureg.convert(1, 'cP', 'Pa*s')
+                * m.dynamic_viscosity
                 * ureg.convert(m.vel_ub[sol], 'cm', 'm')
                 / m.dh[sol] ** 2
                 * m.flow_vol['rm', sol].ub
@@ -1688,7 +1694,7 @@ def _deltaP(m, x, sol):
         """
         return ureg.convert(
             48
-            * ureg.convert(1, 'cP', 'Pa*s')
+            * m.dynamic_viscosity
             * ureg.convert(m.vel[x, sol], 'cm', 'm')
             / m.dh[sol] ** 2,
             'Pa',