From c7665bc79c5ca78c430198b412af19ef1be4f96b Mon Sep 17 00:00:00 2001 From: milanofthe Date: Thu, 7 May 2026 10:57:29 +0200 Subject: [PATCH 1/4] =?UTF-8?q?Fix=20flash=20drum=20holdup=20dynamics=20?= =?UTF-8?q?=E2=80=94=20VLE=20on=20drum=20liquid,=20not=20feed?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/pathsim_chem/process/flash_drum.py | 34 ++++++++--- .../process/multicomponent_flash.py | 34 ++++++++--- tests/process/test_flash_drum.py | 56 +++++++++++++++++++ 3 files changed, 106 insertions(+), 18 deletions(-) diff --git a/src/pathsim_chem/process/flash_drum.py b/src/pathsim_chem/process/flash_drum.py index 29f71a0..0c94513 100644 --- a/src/pathsim_chem/process/flash_drum.py +++ b/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, diff --git a/src/pathsim_chem/process/multicomponent_flash.py b/src/pathsim_chem/process/multicomponent_flash.py index a5e3b24..4cbbe2a 100644 --- a/src/pathsim_chem/process/multicomponent_flash.py +++ b/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,7 +216,7 @@ 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] @@ -215,14 +224,18 @@ def _fn_d(x, u, t): 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 diff --git a/tests/process/test_flash_drum.py b/tests/process/test_flash_drum.py index 454e04a..37fe934 100644 --- a/tests/process/test_flash_drum.py +++ b/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 ==================================================================== From 0cbfa12d7c799244cd16ab6767686f4d16e68e4c Mon Sep 17 00:00:00 2001 From: milanofthe Date: Thu, 7 May 2026 10:58:46 +0200 Subject: [PATCH 2/4] Fix UNIQUAC NaN at pure component using algebraic V/x and F/V identities --- .../thermodynamics/activity_coefficients.py | 15 ++++++++++----- .../thermodynamics/test_activity_coefficients.py | 3 +++ 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/pathsim_chem/thermodynamics/activity_coefficients.py b/src/pathsim_chem/thermodynamics/activity_coefficients.py index a4ca198..bd73821 100644 --- a/src/pathsim_chem/thermodynamics/activity_coefficients.py +++ b/src/pathsim_chem/thermodynamics/activity_coefficients.py @@ -306,17 +306,22 @@ def _eval(self, T): tau = np.exp(self.a + self.b_param / T + self.c_param * np.log(T) + self.d_param * T) # volume and surface fractions - V = self.r * x / np.dot(self.r, x) - F = self.q * x / np.dot(self.q, x) + sum_rx = np.dot(self.r, x) + sum_qx = np.dot(self.q, x) Fp = self.q_prime * x / np.dot(self.q_prime, x) # combinatorial part + # V_i / x_i = r_i / sum_j(r_j x_j) is well-defined even at x_i=0 + # F_i / V_i = (q_i / r_i) * sum_rx / sum_qx (algebraic identity) + V_over_x = self.r / sum_rx + F_over_V = (self.q / self.r) * (sum_rx / sum_qx) + ln_gamma_C = np.zeros(n) sum_xl = np.dot(x, self.l) for i in range(n): - ln_gamma_C[i] = (np.log(V[i] / x[i]) - + self.z / 2 * self.q[i] * np.log(F[i] / V[i]) - + self.l[i] - V[i] / x[i] * sum_xl) + ln_gamma_C[i] = (np.log(V_over_x[i]) + + self.z / 2 * self.q[i] * np.log(F_over_V[i]) + + self.l[i] - V_over_x[i] * sum_xl) # residual part ln_gamma_R = np.zeros(n) diff --git a/tests/thermodynamics/test_activity_coefficients.py b/tests/thermodynamics/test_activity_coefficients.py index f6ad11b..aec047d 100644 --- a/tests/thermodynamics/test_activity_coefficients.py +++ b/tests/thermodynamics/test_activity_coefficients.py @@ -174,6 +174,9 @@ def test_pure_component_gamma_is_one(self): ) gamma1, gamma2 = eval_block(U, 350) self.assertAlmostEqual(gamma1, 1.0, places=10) + # gamma_2 at infinite dilution must be finite (not NaN) + self.assertTrue(np.isfinite(gamma2)) + self.assertGreater(gamma2, 0.0) def test_symmetric_case(self): # Identical r, q, symmetric a => gamma_1 = gamma_2 From 16dea377bee8ebf76e2cb33238627d519858959e Mon Sep 17 00:00:00 2001 From: milanofthe Date: Thu, 7 May 2026 11:00:09 +0200 Subject: [PATCH 3/4] Fix Redlich-Kister excess enthalpy: standard polynomial form, drop spurious 0.5 factor --- src/pathsim_chem/thermodynamics/enthalpy.py | 31 +++++++++++--------- tests/thermodynamics/test_enthalpy.py | 32 ++++++++++++++------- 2 files changed, 40 insertions(+), 23 deletions(-) diff --git a/src/pathsim_chem/thermodynamics/enthalpy.py b/src/pathsim_chem/thermodynamics/enthalpy.py index 24ff03e..f8de53f 100644 --- a/src/pathsim_chem/thermodynamics/enthalpy.py +++ b/src/pathsim_chem/thermodynamics/enthalpy.py @@ -301,22 +301,24 @@ class ExcessEnthalpyRedlichKister(Function): Computes the molar excess enthalpy using a Redlich-Kister polynomial expansion. This is a flexible, empirical model for representing binary excess properties. For a multi-component mixture, the excess enthalpy - is computed as a sum of binary pair contributions: + is summed over the unique binary pairs supplied: .. math:: - h^E = \frac{1}{2} \sum_i \sum_j h^E_{ij} + h^E = \sum_{(i,j)} h^E_{ij} - Each binary pair contribution: + Each binary pair contribution follows the standard Redlich-Kister form .. math:: h^E_{ij} = \frac{x_i x_j}{x_i + x_j} - \left(A(T) x_d + B(T) x_d^2 + C(T) x_d^3 + \ldots\right) + \left(A_0(T) + A_1(T) x_d + A_2(T) x_d^2 + \ldots\right) - where :math:`x_d = x_i - x_j` and the coefficients :math:`A(T)`, - :math:`B(T)`, ... are temperature-dependent polynomials: - :math:`A(T) = a_0 + a_1 T + a_2 T^2 + \ldots` + where :math:`x_d = x_i - x_j` and each :math:`A_k(T)` is a temperature + polynomial :math:`A_k(T) = a_{k,0} + a_{k,1} T + a_{k,2} T^2 + \ldots`. + For a true binary system :math:`x_i + x_j = 1` and the prefactor + reduces to :math:`x_i x_j`; in multicomponent mixtures the + :math:`/(x_i+x_j)` term provides the symmetric extension. **Input port:** ``T`` -- temperature [K]. @@ -328,12 +330,15 @@ class ExcessEnthalpyRedlichKister(Function): Mole fractions [N]. coeffs : dict Redlich-Kister coefficients keyed by binary pair ``(i, j)`` as tuples. - Each value is a list of polynomial coefficient arrays, one per - Redlich-Kister term. Example for a single pair (0,1) with 3 terms:: + Provide each unordered pair only once (e.g. ``(0, 1)``, not also + ``(1, 0)``). Each value is a list of polynomial coefficient arrays, + one per Redlich-Kister term. Example for a single pair (0,1) with + 3 terms:: {(0, 1): [[a0, a1], [b0, b1], [c0]]} - This gives A(T)=a0+a1*T, B(T)=b0+b1*T, C(T)=c0. + This gives A_0(T)=a0+a1*T, A_1(T)=b0+b1*T, A_2(T)=c0, where + :math:`A_k` multiplies :math:`x_d^k`. """ input_port_labels = {"T": 0} @@ -358,9 +363,9 @@ def _eval(self, T): xd = xi - xj - # evaluate each Redlich-Kister term + # evaluate Redlich-Kister polynomial: A_0 + A_1*xd + A_2*xd^2 + ... pair_hE = 0.0 - xd_power = xd + xd_power = 1.0 for poly_coeffs in terms: # temperature polynomial: c0 + c1*T + c2*T^2 + ... coeff_val = 0.0 @@ -373,7 +378,7 @@ def _eval(self, T): hE += xi * xj / x_sum * pair_hE - return 0.5 * hE + return hE # ENTHALPY DEPARTURE FROM EOS (9.7) ===================================================== diff --git a/tests/thermodynamics/test_enthalpy.py b/tests/thermodynamics/test_enthalpy.py index 09c8740..a1f0dd2 100644 --- a/tests/thermodynamics/test_enthalpy.py +++ b/tests/thermodynamics/test_enthalpy.py @@ -201,28 +201,40 @@ def test_init(self): ) self.assertEqual(hE.n, 2) - def test_symmetric_binary(self): - # Simplest case: one-term Redlich-Kister with constant A - # h^E_ij = x_i*x_j/(x_i+x_j) * A * (x_i-x_j) - # For x=[0.5, 0.5]: h^E = 0.5 * 0.5 / 1.0 * A * 0 = 0 + def test_constant_coeff_only(self): + # One-term Redlich-Kister with constant A_0: + # h^E = x_i*x_j/(x_i+x_j) * A_0 + # For x=[0.5, 0.5], A_0=1000: h^E = 0.25/1.0 * 1000 = 250 hE = ExcessEnthalpyRedlichKister( x=[0.5, 0.5], coeffs={(0, 1): [[1000.0]]}, ) result = eval_block_T(hE, 300) - self.assertAlmostEqual(result, 0.0, places=10) + self.assertAlmostEqual(result, 250.0, places=10) - def test_asymmetric_composition(self): - # x=[0.3, 0.7], one-term A=2000 - # h^E = 0.5 * 0.3*0.7/1.0 * 2000 * (0.3-0.7) = 0.5 * 0.21 * 2000 * (-0.4) = -84 + def test_two_term_polynomial(self): + # x=[0.3, 0.7], coeffs A_0=2000, A_1=500 + # h^E = 0.21 * (2000 + 500*(-0.4)) = 0.21 * 1800 = 378 hE = ExcessEnthalpyRedlichKister( x=[0.3, 0.7], - coeffs={(0, 1): [[2000.0]]}, + coeffs={(0, 1): [[2000.0], [500.0]]}, ) result = eval_block_T(hE, 300) - expected = 0.5 * 0.3 * 0.7 / 1.0 * 2000.0 * (0.3 - 0.7) + expected = 0.3 * 0.7 / 1.0 * (2000.0 + 500.0 * (0.3 - 0.7)) self.assertAlmostEqual(result, expected, places=5) + def test_antisymmetric_under_pair_swap(self): + # Odd-order RK terms should flip sign when (i,j) <-> (j,i), + # since x_d = x_i - x_j flips. Verify by swapping x values. + hE_a = ExcessEnthalpyRedlichKister( + x=[0.3, 0.7], coeffs={(0, 1): [[0.0], [1000.0]]}, # A_1 only + ) + hE_b = ExcessEnthalpyRedlichKister( + x=[0.7, 0.3], coeffs={(0, 1): [[0.0], [1000.0]]}, + ) + self.assertAlmostEqual(eval_block_T(hE_a, 300), + -eval_block_T(hE_b, 300), places=8) + def test_temperature_dependent_coeffs(self): # A(T) = 1000 + 2*T hE = ExcessEnthalpyRedlichKister( From 104f9203cbe6378f071d3736edfb7ee6acf95033 Mon Sep 17 00:00:00 2001 From: milanofthe Date: Thu, 7 May 2026 11:00:31 +0200 Subject: [PATCH 4/4] Move TCAP1D NotImplementedError to __init__ so import succeeds --- src/pathsim_chem/tritium/tcap.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pathsim_chem/tritium/tcap.py b/src/pathsim_chem/tritium/tcap.py index 06df767..2825594 100644 --- a/src/pathsim_chem/tritium/tcap.py +++ b/src/pathsim_chem/tritium/tcap.py @@ -12,15 +12,16 @@ # BLOCKS ================================================================================ class TCAP1D(ODE): - """This block models the Thermal Cycle Absorption Process (TCAP) in 1d. + """This block models the Thermal Cycle Absorption Process (TCAP) in 1d. - The model uses a 1d finite difference spatial discretization to construct + The model uses a 1d finite difference spatial discretization to construct a nonlinear ODE internally as proposed in https://doi.org/10.1016/j.ijhydene.2023.03.101 - """ - raise NotImplementedError("TCAP1D block is currently not impolemented!") + + def __init__(self, *args, **kwargs): + raise NotImplementedError("TCAP1D block is not yet implemented")