Fix flash drum holdup dynamics — VLE on drum liquid, not feed

Author: milanofthe

src/pathsim_chem/process/flash_drum.py

@@ -43,11 +43,20 @@ class FlashDrum(DynamicalSystem):
 
     Dynamic States
     ---------------
-    The holdup moles of each component in the liquid phase:
+    The drum holds well-mixed liquid; vapor exits immediately at feed-flash
+    composition. Liquid component holdup follows a CSTR balance with the
+    feed-flash liquid as inlet:
 
     .. math::
 
-        \\frac{dN_i}{dt} = F z_i - V y_i - L x_i
+        \\frac{dN_i}{dt} = L_{in} \\, (x_{eq,i} - x_{drum,i})
+
+    where :math:`L_{in} = (1-\\beta) F` is the liquid feed to the drum,
+    :math:`x_{eq}` the Rachford-Rice liquid composition for the feed, and
+    :math:`x_{drum} = N / \\sum_j N_j` the drum liquid composition. Total
+    holdup :math:`M = \\sum_i N_i` is preserved exactly. The output ``x_1``
+    is :math:`x_{drum,1}` (dynamic), while ``y_1``, ``V_rate``, ``L_rate``
+    follow the instantaneous feed flash.
 
     Parameters
     ----------
@@ -147,20 +156,24 @@ def _pad_u(u):
                 return padded
             return u
 
-        # rhs of flash drum ode
+        # rhs of flash drum ode: liquid CSTR fed by feed-flash liquid
         def _fn_d(x, u, t):
             u = _pad_u(u)
             F_in, z_1, T, P = u
             z = np.array([z_1, 1.0 - z_1])
 
-            beta, y, x_eq = _solve_vle(z, T, P)
+            if F_in <= 0.0:
+                return np.zeros(2)
 
-            V_rate = beta * F_in
-            L_rate = (1.0 - beta) * F_in
+            beta, _, x_eq = _solve_vle(z, T, P)
+            L_in = (1.0 - beta) * F_in
+
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
 
-            return F_in * z - V_rate * y - L_rate * x_eq
+            return L_in * (x_eq - x_drum)
 
-        # algebraic output
+        # algebraic output: V/L/y from feed flash (instantaneous), x from drum
         def _fn_a(x, u, t):
             u = _pad_u(u)
             F_in, z_1, T, P = u
@@ -171,7 +184,10 @@ def _fn_a(x, u, t):
             V_rate = beta * F_in
             L_rate = (1.0 - beta) * F_in
 
-            return np.array([V_rate, L_rate, y[0], x_eq[0]])
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
+
+            return np.array([V_rate, L_rate, y[0], x_drum[0]])
 
         super().__init__(
             func_dyn=_fn_d,

src/pathsim_chem/process/multicomponent_flash.py

@@ -35,11 +35,20 @@ class MultiComponentFlash(DynamicalSystem):
 
     Dynamic States
     ---------------
-    The holdup moles of each component in the liquid phase:
+    The drum holds well-mixed liquid; vapor exits immediately at feed-flash
+    composition. Liquid component holdup follows a CSTR balance with the
+    feed-flash liquid as inlet:
 
     .. math::
 
-        \\frac{dN_i}{dt} = F z_i - V y_i - L x_i
+        \\frac{dN_i}{dt} = L_{in} \\, (x_{eq,i} - x_{drum,i})
+
+    where :math:`L_{in} = (1-\\beta) F` is the liquid feed to the drum,
+    :math:`x_{eq}` the Rachford-Rice liquid composition for the feed, and
+    :math:`x_{drum} = N / \\sum_j N_j` the drum liquid composition. Total
+    holdup :math:`M = \\sum_i N_i` is preserved exactly. Outputs ``x_i``
+    are :math:`x_{drum,i}` (dynamic); ``V_rate``, ``L_rate``, ``y_i`` come
+    from the instantaneous feed flash.
 
     Parameters
     ----------
@@ -207,22 +216,26 @@ def _extract_z(u):
             z_last = 1.0 - np.sum(z_partial)
             return np.append(z_partial, z_last)
 
-        # rhs of flash drum ode
+        # rhs of flash drum ode: liquid CSTR fed by feed-flash liquid
         def _fn_d(x, u, t):
             u = _pad_u(u)
             F_in = u[0]
             z = _extract_z(u)
             T = u[nc]
             P = u[nc + 1]
 
-            beta, y, x_eq = _solve_vle(z, T, P)
+            if F_in <= 0.0:
+                return np.zeros(nc)
 
-            V_rate = beta * F_in
-            L_rate = (1.0 - beta) * F_in
+            beta, _, x_eq = _solve_vle(z, T, P)
+            L_in = (1.0 - beta) * F_in
+
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
 
-            return F_in * z - V_rate * y - L_rate * x_eq
+            return L_in * (x_eq - x_drum)
 
-        # algebraic output
+        # algebraic output: V/L/y from feed flash (instant), x from drum state
         def _fn_a(x, u, t):
             u = _pad_u(u)
             F_in = u[0]
@@ -235,12 +248,15 @@ def _fn_a(x, u, t):
             V_rate = beta * F_in
             L_rate = (1.0 - beta) * F_in
 
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
+
             # output: V_rate, L_rate, y_1..y_{nc-1}, x_1..x_{nc-1}
             result = np.empty(2 + 2 * (nc - 1))
             result[0] = V_rate
             result[1] = L_rate
             result[2:2 + nc - 1] = y[:nc - 1]
-            result[2 + nc - 1:] = x_eq[:nc - 1]
+            result[2 + nc - 1:] = x_drum[:nc - 1]
 
             return result
 

tests/process/test_flash_drum.py

@@ -136,6 +136,62 @@ def test_no_feed(self):
         self.assertAlmostEqual(F.outputs[0], 0.0)  # V_rate
         self.assertAlmostEqual(F.outputs[1], 0.0)  # L_rate
 
+    def test_holdup_dynamics_drives_to_equilibrium(self):
+        """At steady state the drum liquid composition must equal the RR
+        equilibrium liquid composition for the feed."""
+        F = FlashDrum(holdup=100.0, N0=[80.0, 20.0])  # off-equilibrium init
+        F.set_solver(EUF, parent=None)
+
+        # T=370 K gives two-phase region for benzene/toluene defaults at 1 atm
+        T, P = 370.0, 101325.0
+        u = np.array([10.0, 0.5, T, P])
+
+        # at the RR equilibrium x_eq, dN/dt must vanish
+        # compute x_eq via direct VLE (binary RR with same Antoine defaults)
+        Psat = np.exp(F.antoine_A - F.antoine_B / (T + F.antoine_C))
+        K = Psat / P
+        z = np.array([0.5, 0.5])
+        d1, d2 = K[0] - 1, K[1] - 1
+        beta = -(z[0]*d1 + z[1]*d2) / (d1*d2)
+        x_eq = z / (1.0 + beta * (K - 1.0))
+        x_eq = x_eq / x_eq.sum()
+
+        # state at equilibrium with same total holdup
+        N_eq = 100.0 * x_eq
+        dN = F.op_dyn(N_eq, u, 0.0)
+        self.assertTrue(np.allclose(dN, 0.0, atol=1e-10))
+
+        # state away from equilibrium: dN must be non-zero
+        dN_off = F.op_dyn(np.array([80.0, 20.0]), u, 0.0)
+        self.assertGreater(np.linalg.norm(dN_off), 1e-3)
+
+    def test_holdup_total_moles_conserved(self):
+        """dM/dt = sum(dN/dt) must be exactly zero (perfect level control)."""
+        F = FlashDrum(holdup=100.0, N0=[70.0, 30.0])
+        F.set_solver(EUF, parent=None)
+
+        u = np.array([5.0, 0.4, 355.0, 101325.0])
+        for state in (np.array([70.0, 30.0]),
+                      np.array([20.0, 80.0]),
+                      np.array([1.0, 99.0])):
+            dN = F.op_dyn(state, u, 0.0)
+            self.assertAlmostEqual(dN.sum(), 0.0, places=10,
+                                    msg=f"dM/dt != 0 for state {state}")
+
+    def test_x_output_uses_drum_state(self):
+        """x_1 output must reflect drum state, not feed composition."""
+        F = FlashDrum(holdup=100.0, N0=[90.0, 10.0])
+        F.set_solver(EUF, parent=None)
+
+        F.inputs[0] = 10.0
+        F.inputs[1] = 0.3       # different from drum
+        F.inputs[2] = 360.0
+        F.inputs[3] = 101325.0
+
+        F.update(None)
+        # drum state x_1 = 90/100 = 0.9
+        self.assertAlmostEqual(F.outputs[3], 0.9, places=8)
+
 
 # RUN TESTS LOCALLY ====================================================================