Add implementation for simple dead end filtation unit

Author: schmoelder

CADETPythonSimulator/residual.py

@@ -77,6 +77,70 @@ def calculate_residual_visc_cstr():
     return 0
 
 
-def calculate_residual_def():
-    """Calculate the residual equations fo a dead end filtration equation."""
-    raise NotImplementedError
+def calculate_residual_cake_vol_def(
+        V_dot_f: float,
+        rejection: np.ndarray,
+        molar_volume: np.ndarray,
+        c_in: np.ndarray,
+        V_dot_C: float
+        ) -> float:
+    """
+    Residual equation for the Cake Volume.
+
+    Parameters
+    ----------
+    V_dot_f : float
+        Flowrate of incoming feed
+    rejection : float
+        Rejection of the filter
+    gamma : float
+        Portion of suspended material
+    V_dot_C : float
+        Change of Cake Volume
+    """
+    return -V_dot_C + np.sum(rejection * molar_volume * c_in * V_dot_f)
+
+
+def calculate_residual_press_easy_def(
+        V_dot_Perm: float,
+        V_C: float,
+        deltap: float,
+        A: float,
+        mu: float,
+        Rm: float,
+        alpha: float
+        ) -> float:
+    """
+    Calculate the residual equations fo a dead end filtration equation for the pressure
+    in the easy model.
+
+    Parameters
+    ----------
+    V_dot_Perm : np.ndarray
+        FLow of the Permeate through the membrane and Cake
+    V_C : float
+        Volume of the Cake
+    deltap : float
+        Pressure drop in this unit
+    A : float
+        Filtration area
+    mu : float
+        dynamic Viscosity
+    Rm : float
+        resistance of the medium
+    alpha : float
+        Specific cake resistance
+    """
+    hyd_resistance = (Rm + alpha*V_C/A) * mu
+
+    return -V_dot_Perm + deltap * A *hyd_resistance
+
+
+
+def calculate_residual_visc_def():
+    """
+    Calculate the residual of the Viscosity equation of the CSTR.
+    """
+    warnings.warn("Viscosity of def not yet implemented")
+
+    return 0

CADETPythonSimulator/unit_operation.py

@@ -15,7 +15,12 @@
 from CADETPythonSimulator.exception import NotInitializedError, CADETPythonSimError
 from CADETPythonSimulator.state import State, state_factory
 from CADETPythonSimulator.residual import (
-    calculate_residual_volume_cstr, calculate_residual_concentration_cstr, calculate_residuals_visc_cstr
+    calculate_residual_volume_cstr,
+    calculate_residual_concentration_cstr,
+    calculate_residual_visc_cstr,
+    calculate_residual_press_easy_def,
+    calculate_residual_cake_vol_def,
+    calculate_residual_visc_def
     )
 from CADETPythonSimulator.rejection import RejectionBase
 from CADETPythonSimulator.cake_compressibility import CakeCompressibilityBase
@@ -593,73 +598,121 @@ class DeadEndFiltration(UnitOperationBase):
         Model for cake compressibility.
 
     """
-    retentate = {
+    cake = {
         'dimensions': (),
-        'entries': {'c': 'n_comp', 'viscosity': 1, 'Rc': 1, 'mc': 'n_comp'},
+        'entries': {'c': 'n_comp',
+                    'viscosity': 1,
+                    'pressure': 1,
+                    'cakevolume': 1,
+                    'permeate': 1
+                    },
         'n_inlet_ports': 1,
     }
     permeate = {
         'dimensions': (),
         'entries': {'c': 'n_comp', 'viscosity': 1, 'Volume': 1},
         'n_outlet_ports': 1,
     }
-    _state_structures = ['retentate', 'permeate']
-
-    rejection_model = Typed(ty=RejectionBase)
-    cake_compressibility_model = Typed(ty=CakeCompressibilityBase)
+    _state_structures = ['cake', 'permeate']
 
     membrane_area = UnsignedFloat()
     membrane_resistance = UnsignedFloat()
+    specific_cake_resistance = UnsignedFloat()
+    rejection = Typed(ty=RejectionBase)
 
     _parameters = [
         'membrane_area',
         'membrane_resistance',
+        'specific_cake_resistance',
+        'rejection'
     ]
 
-    def delta_p(self):
-        raise NotImplementedError()
+    def compute_residual(
+            self,
+            t: float,
+            ) -> NoReturn:
 
-    def specific_cake_resistance(self, delta_p: float) -> float:
-        """
-        Compute specific resistance as a function of delta_p.
+        Q_in = self.Q_in[0]
+        Q_out = self.Q_out[0]
 
-        Parameters
-        ----------
-        delta_p : float
-            Pressure difference.
+        c_in = self.states['cake']['c']
+        c_in_dot = self.state_derivatives['cake']['c']
 
-        Returns
-        -------
-        float
-            Specific cake resistance.
+        V_C = self.states['cake']['cakevolume']
+        V_dot_C = self.state_derivatives['cake']['cakevolume']
 
-        """
-        raise self.cake_compressibility_model.specific_cake_resistance(delta_p)
+        V_p = self.states['cake']['permeate']
+        Q_p = self.state_derivatives['cake']['cakevolume']
 
-    def compute_residual(
-            self,
-            t: float,
-            y: np.ndarray,
-            y_dot: np.ndarray,
-            residual: np.ndarray
-            ) -> NoReturn:
-        #     0,  1,  2
-        # y = Vp, Rc, mc
-        # TODO: Needs to be extended to include c_in / c_out
-        # y = [*c_i_in], viscosity_in, Vp, Rc, mc, [*c_i_out], viscosity_out
+        viscosity_in = self.states['cake']['viscosity']
+
+        c = self.states['permeate']['c']
+        c_dot = self.state_derivatives['permeate']['c']
 
-        c_in = y[0: self.n_comp]
-        viscosity_in = y[self.n_comp]
+        V = self.states['permeate']['Volume']
+        V_dot = self.state_derivatives['permeate']['Volume']
 
-        densities = self.component_system.densities
+        deltap = self.states['cake']['pressure']
 
-        residual[self.n_dof_coupling + 0] = ((self.membrane_area*self.delta_p(t)/viscosity_in)/(self.membrane_resistance+y[1])) - y_dot[0]
-        residual[self.n_dof_coupling + 1] = (1/self.membrane_area) * (y_dot[2] * self.specific_cake_resistance(self.p(t))) - y_dot[1]
+        #parameters
+        molecular_weights = self.component_system.molecular_weights
+        molar_volume = self.component_system.molecular_volumes
+        membrane_area = self.parameters['membrane_area']
+        membrane_resistance = self.parameters['membrane_resistance']
+        specific_cake_resistance = self.parameters['specific_cake_resistance']
 
-        residual[self.n_dof_coupling + 2] = ((self.c(t) * y_dot[0]) / (1-self.c(t)/self.density)) - y_dot[2]
+        rejection = np.array(
+                        [self.rejection.get_rejection(mw) for mw in molecular_weights])
 
-        self.residuals['retentate']
-        self.residuals['permeate']
+        # Handle inlet DOFs, which are simply copied to the residual
+        self.residuals['cake']['c'] = c_in
+        self.residuals['cake']['cakevolume'] = calculate_residual_cake_vol_def(
+                                                Q_in,
+                                                rejection,
+                                                molar_volume,
+                                                c_in,
+                                                V_dot_C
+                                                )
+
+        self.residuals['cake']['pressure'] = calculate_residual_press_easy_def(
+                                                Q_p,
+                                                V_C,
+                                                deltap,
+                                                membrane_area,
+                                                viscosity_in,
+                                                membrane_resistance,
+                                                specific_cake_resistance
+                                                )
+
+        self.residuals['cake']['permeate'] = calculate_residual_volume_cstr(
+                                                V_C,
+                                                V_dot_C,
+                                                Q_in,
+                                                Q_p
+                                                )
+
+        self.residuals['cake']['viscosity'] = calculate_residual_visc_def()
+
+        new_c_in = (1-rejection)*c_in
+
+        self.residuals['permeate']['c'] = calculate_residual_concentration_cstr(
+                                                c,
+                                                c_dot,
+                                                V,
+                                                V_dot,
+                                                Q_p,
+                                                Q_out,
+                                                new_c_in
+                                                )
+
+        self.residuals['permeate']['Volume'] = calculate_residual_volume_cstr(
+                                                V,
+                                                V_dot,
+                                                Q_p,
+                                                Q_out
+                                                )
+
+        self.residuals['permeate']['viscosity'] = calculate_residual_visc_cstr()
 
 
 

tests/test_residual.py

@@ -2,7 +2,10 @@
 import pytest
 
 from CADETPythonSimulator.residual import (
-    calculate_residual_volume_cstr, calculate_residual_concentration_cstr
+    calculate_residual_volume_cstr,
+    calculate_residual_concentration_cstr,
+    calculate_residual_cake_vol_def,
+    calculate_residual_press_easy_def
 )
 from CADETPythonSimulator.exception import CADETPythonSimError
 
@@ -172,6 +175,114 @@ def test_calculation_cstr(self, parameters):
 
         residual = calculate_residual_volume_cstr(*param_vec_volume)
 
+# Testcase 1: Membrane rejects all
+TestCaseDEFCake_rejects_all = {
+    "values": {
+        "V_dot_f": 1.0,
+        "rejection": np.array([1, 1]),
+        "molar_volume": np.array([1, 1]),
+        "c_in": np.array([0.5, 0.5]),
+        "V_dot_C": 1.0
+    },
+    "expected": 0
+}
+
+
+# Testcase 2: Membrane rejects nothing
+TestCaseDEFCake_rejects_not = {
+    "values": {
+        "V_dot_f": 1.0,
+        "rejection": np.array([0, 0]),
+        "molar_volume": np.array([1, 1]),
+        "c_in": np.array([0.5, 0.5]),
+        "V_dot_C": 0.0
+    },
+    "expected": 0
+}
+
+# Testcase 3: Membrane rejects only Component 2
+TestCaseDEFCake_rejects_2 = {
+    "values": {
+        "V_dot_f": 1.0,
+        "rejection": np.array([0, 1]),
+        "molar_volume": np.array([1, 1]),
+        "c_in": np.array([0.5, 0.5]),
+        "V_dot_C": 0.5
+    },
+    "expected": 0
+}
+
+# Testcase 4: Component 2 is larger then 1
+TestCaseDEFCake_C2_le_C1 = {
+    "values": {
+        "V_dot_f": 1.0,
+        "rejection": np.array([1, 1]),
+        "molar_volume": np.array([0.5, 1]),
+        "c_in": np.array([0.5, 0.5]),
+        "V_dot_C": 0.75
+    },
+    "expected": 0
+}
+
+
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        TestCaseDEFCake_rejects_all,
+        TestCaseDEFCake_rejects_not,
+        TestCaseDEFCake_rejects_2,
+        TestCaseDEFCake_C2_le_C1
+    ]
+)
+class TestResidualCakeVolDEF():
+    def test_calculation_def(self, parameters):
+        param_vec_cake_vol = parameters["values"].values()
+        np.testing.assert_equal(
+            calculate_residual_cake_vol_def(*param_vec_cake_vol),
+            parameters["expected"]
+        )
+
+
+# Case 1 : Equally large hyraulic resistance
+TestCaseDEFPressureDrop = {
+    "values": {
+        "V_dot_P": 1,
+        "V_C": 1,
+        "deltap": 0.5,
+        "A": 1,
+        "mu": 1,
+        "Rm": 1,
+        "alpha": 1,
+    },
+    "expected": 0
+}
+
+# Case 2 : No cake yet
+TestCaseDEFPressureDrop_no_cake = {
+    "values": {
+        "V_dot_P": 0.5,
+        "V_C": 0,
+        "deltap": 0.5,
+        "A": 1,
+        "mu": 1,
+        "Rm": 1,
+        "alpha": 1,
+    },
+    "expected": 0
+}
+
+
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        TestCaseDEFPressureDrop,
+        TestCaseDEFPressureDrop_no_cake
+    ]
+)
+class TestResidualPressureDropDEF():
+    def test_calculation_def(self, parameters):
+        param_vec_pressure = parameters["values"].values()
+        residual = calculate_residual_press_easy_def(*param_vec_pressure)
         np.testing.assert_equal(residual, parameters["expected"])
 
 
@@ -198,14 +309,12 @@ def test_calculation_cstr(self, parameters):
 class TestResidualError():
 
     def test_calculation_vol_cstr_error(self, parameters):
-
         param_vec_volume = parameters["values"].values()
 
         with pytest.raises(CADETPythonSimError):
             calculate_residual_volume_cstr(*list(param_vec_volume)[2:6])
 
     def test_calculation_concentration_cstr_error(self, parameters):
-
         param_vec_volume = parameters["values"].values()
 
         with pytest.raises(CADETPythonSimError):