From 034015cf156ed101ca3b428e7d73a5b29001fed7 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 08:44:03 +0100 Subject: [PATCH 1/9] Add IK-CAPE pure component property correlations (16 blocks) --- src/pathsim_chem/__init__.py | 1 + src/pathsim_chem/thermodynamics/__init__.py | 1 + .../thermodynamics/correlations.py | 477 ++++++++++++++++++ tests/thermodynamics/__init__.py | 0 tests/thermodynamics/test_correlations.py | 383 ++++++++++++++ 5 files changed, 862 insertions(+) create mode 100644 src/pathsim_chem/thermodynamics/__init__.py create mode 100644 src/pathsim_chem/thermodynamics/correlations.py create mode 100644 tests/thermodynamics/__init__.py create mode 100644 tests/thermodynamics/test_correlations.py diff --git a/src/pathsim_chem/__init__.py b/src/pathsim_chem/__init__.py index 1f4f7c6..a879933 100644 --- a/src/pathsim_chem/__init__.py +++ b/src/pathsim_chem/__init__.py @@ -14,3 +14,4 @@ #for direct block import from main package from .tritium import * +from .thermodynamics import * diff --git a/src/pathsim_chem/thermodynamics/__init__.py b/src/pathsim_chem/thermodynamics/__init__.py new file mode 100644 index 0000000..f4a6307 --- /dev/null +++ b/src/pathsim_chem/thermodynamics/__init__.py @@ -0,0 +1 @@ +from .correlations import * diff --git a/src/pathsim_chem/thermodynamics/correlations.py b/src/pathsim_chem/thermodynamics/correlations.py new file mode 100644 index 0000000..96f25fe --- /dev/null +++ b/src/pathsim_chem/thermodynamics/correlations.py @@ -0,0 +1,477 @@ +######################################################################################### +## +## IK-CAPE Pure Component Property Correlations +## +## Temperature-dependent correlation functions for thermodynamic properties +## as defined in the IK-CAPE Thermodynamics Module, Chapter 2. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# BLOCKS ================================================================================ + +class Polynomial(Function): + r"""Polynomial correlation (IK-CAPE code: POLY). + + .. math:: + + f(T) = \sum_{i=0}^{n} a_i \, T^i + + General-purpose polynomial for temperature-dependent properties. + + Parameters + ---------- + coeffs : list or tuple of float + Polynomial coefficients :math:`[a_0, a_1, \ldots, a_n]` (up to 10). + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, coeffs): + self.coeffs = tuple(coeffs) + super().__init__(func=self._eval) + + def _eval(self, T): + return sum(a * T**i for i, a in enumerate(self.coeffs)) + + +class ExponentialPolynomial(Function): + r"""Exponential polynomial correlation (IK-CAPE code: EPOL). + + .. math:: + + f(T) = 10^{\sum_{i=0}^{n} a_i \, T^i} + + General-purpose correlation for properties that vary over orders of magnitude. + + Parameters + ---------- + coeffs : list or tuple of float + Polynomial exponent coefficients :math:`[a_0, a_1, \ldots, a_n]` (up to 10). + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, coeffs): + self.coeffs = tuple(coeffs) + super().__init__(func=self._eval) + + def _eval(self, T): + exponent = sum(a * T**i for i, a in enumerate(self.coeffs)) + return 10**exponent + + +class Watson(Function): + r"""Watson correlation (IK-CAPE code: WATS). + + .. math:: + + f(T) = a_0 \, (a_2 - T)^{a_1} + a_3 + + Commonly used for heat of vaporization. + + Parameters + ---------- + a0, a1, a2, a3 : float + Watson correlation coefficients. :math:`a_2` is typically the + critical temperature. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2, a3=0.0): + self.coeffs = (a0, a1, a2, a3) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3 = self.coeffs + return a0 * (a2 - T)**a1 + a3 + + +class Antoine(Function): + r"""Antoine correlation (IK-CAPE code: ANTO). + + .. math:: + + \ln f(T) = a_0 - \frac{a_1}{T + a_2} + + Standard three-parameter vapor pressure correlation. + + Parameters + ---------- + a0, a1, a2 : float + Antoine coefficients. For natural-log form with T in Kelvin. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2): + self.coeffs = (a0, a1, a2) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2 = self.coeffs + return np.exp(a0 - a1 / (T + a2)) + + +class ExtendedAntoine(Function): + r"""Extended Antoine correlation (IK-CAPE code: ANT1). + + .. math:: + + \ln f(T) = a_0 + \frac{a_1}{T + a_2} + a_3 \, T + a_4 \, \ln(T) + a_5 \, T^{a_6} + + Extended form for wider temperature range vapor pressure. + + Parameters + ---------- + a0, a1, a2, a3, a4, a5, a6 : float + Extended Antoine coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2, a3=0.0, a4=0.0, a5=0.0, a6=0.0): + self.coeffs = (a0, a1, a2, a3, a4, a5, a6) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4, a5, a6 = self.coeffs + ln_f = a0 + a1 / (T + a2) + a3 * T + a4 * np.log(T) + a5 * T**a6 + return np.exp(ln_f) + + +class Kirchhoff(Function): + r"""Kirchhoff correlation (IK-CAPE code: KIRC). + + .. math:: + + \ln f(T) = a_0 - \frac{a_1}{T} - a_2 \, \ln(T) + + Three-parameter vapor pressure correlation. + + Parameters + ---------- + a0, a1, a2 : float + Kirchhoff coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2): + self.coeffs = (a0, a1, a2) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2 = self.coeffs + return np.exp(a0 - a1 / T - a2 * np.log(T)) + + +class ExtendedKirchhoff(Function): + r"""Extended Kirchhoff correlation (IK-CAPE code: KIR1). + + .. math:: + + \ln f(T) = a_0 + \frac{a_1}{T} + a_2 \, \ln(T) + a_3 \, T^{a_4} + + Extended form with additional power term. + + Parameters + ---------- + a0, a1, a2, a3, a4 : float + Extended Kirchhoff coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2, a3=0.0, a4=0.0): + self.coeffs = (a0, a1, a2, a3, a4) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4 = self.coeffs + return np.exp(a0 + a1 / T + a2 * np.log(T) + a3 * T**a4) + + +class Sutherland(Function): + r"""Sutherland correlation (IK-CAPE code: SUTH). + + .. math:: + + f(T) = \frac{a_0 \, \sqrt{T}}{1 + a_1 / T} + + Used for gas-phase viscosity estimation. + + Parameters + ---------- + a0, a1 : float + Sutherland coefficients. :math:`a_1` is the Sutherland constant. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1): + self.coeffs = (a0, a1) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1 = self.coeffs + return a0 * np.sqrt(T) / (1 + a1 / T) + + +class Wagner(Function): + r"""Wagner correlation (IK-CAPE code: WAGN). + + .. math:: + + \ln f(T) = \ln(a_1) + \frac{1}{T_r} \left( a_2 \, \tau + a_3 \, \tau^{1.5} + + a_4 \, \tau^3 + a_5 \, \tau^6 \right) + + where :math:`T_r = T / a_0` and :math:`\tau = 1 - T_r`. + + High-accuracy vapor pressure correlation. + + Parameters + ---------- + Tc : float + Critical temperature :math:`a_0` [K]. + Pc : float + Critical pressure :math:`a_1` [Pa]. + a2, a3, a4, a5 : float + Wagner correlation coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, Tc, Pc, a2, a3, a4, a5): + self.coeffs = (Tc, Pc, a2, a3, a4, a5) + super().__init__(func=self._eval) + + def _eval(self, T): + Tc, Pc, a2, a3, a4, a5 = self.coeffs + Tr = T / Tc + tau = 1 - Tr + ln_f = np.log(Pc) + (1 / Tr) * (a2 * tau + a3 * tau**1.5 + a4 * tau**3 + a5 * tau**6) + return np.exp(ln_f) + + +class LiquidHeatCapacity(Function): + r"""Liquid heat capacity correlation (IK-CAPE code: CPL). + + .. math:: + + f(T) = a_0 + a_1 \, T + a_2 \, T^2 + a_3 \, T^3 + \frac{a_4}{T^2} + + Five-parameter liquid heat capacity correlation. + + Parameters + ---------- + a0, a1, a2, a3, a4 : float + Liquid Cp coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1=0.0, a2=0.0, a3=0.0, a4=0.0): + self.coeffs = (a0, a1, a2, a3, a4) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4 = self.coeffs + return a0 + a1 * T + a2 * T**2 + a3 * T**3 + a4 / T**2 + + +class ExtendedLiquidHeatCapacity(Function): + r"""Extended liquid heat capacity correlation (IK-CAPE code: ICPL). + + .. math:: + + f(T) = a_0 + a_1 \, T + a_2 \, T^2 + a_3 \, T^3 + a_4 \, T^4 + \frac{a_5}{T} + + Six-parameter extended liquid heat capacity correlation. + + Parameters + ---------- + a0, a1, a2, a3, a4, a5 : float + Extended liquid Cp coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1=0.0, a2=0.0, a3=0.0, a4=0.0, a5=0.0): + self.coeffs = (a0, a1, a2, a3, a4, a5) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4, a5 = self.coeffs + return a0 + a1 * T + a2 * T**2 + a3 * T**3 + a4 * T**4 + a5 / T + + +class DynamicViscosity(Function): + r"""Dynamic viscosity correlation (IK-CAPE code: VISC). + + .. math:: + + f(T) = a_0 \, \exp\!\left(\frac{a_1}{T}\right) + a_2 + + Simple two/three-parameter liquid viscosity correlation. + + Parameters + ---------- + a0, a1 : float + Pre-exponential factor and activation energy parameter. + a2 : float, optional + Additive constant (default 0). + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2=0.0): + self.coeffs = (a0, a1, a2) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2 = self.coeffs + return a0 * np.exp(a1 / T) + a2 + + +class Rackett(Function): + r"""Rackett correlation (IK-CAPE code: RACK). + + .. math:: + + f(T) = \frac{a_0}{a_1^{\left(1 + (1 - T/a_2)^{a_3}\right)}} + + Saturated liquid density correlation. + + Parameters + ---------- + a0 : float + Scaling parameter (related to critical volume). + a1 : float + Rackett parameter (base of exponent). + a2 : float + Critical temperature [K]. + a3 : float + Exponent parameter (often 2/7). + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2, a3): + self.coeffs = (a0, a1, a2, a3) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3 = self.coeffs + return a0 / a1**(1 + (1 - T / a2)**a3) + + +class AlyLee(Function): + r"""Aly-Lee correlation (IK-CAPE code: ALYL). + + .. math:: + + f(T) = a_0 + a_1 \left(\frac{a_2 / T}{\sinh(a_2 / T)}\right)^2 + + a_3 \left(\frac{a_4 / T}{\cosh(a_4 / T)}\right)^2 + + Ideal gas heat capacity correlation. + + Parameters + ---------- + a0, a1, a2, a3, a4 : float + Aly-Lee coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2, a3, a4): + self.coeffs = (a0, a1, a2, a3, a4) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4 = self.coeffs + x = a2 / T + y = a4 / T + return a0 + a1 * (x / np.sinh(x))**2 + a3 * (y / np.cosh(y))**2 + + +class DIPPR4(Function): + r"""DIPPR Equation 4 correlation (IK-CAPE code: DIP4). + + .. math:: + + f(T) = a_1 \, (1 - T_r)^h, \quad + h = a_2 + a_3 \, T_r + a_4 \, T_r^2 + a_5 \, T_r^3 + + where :math:`T_r = T / a_0`. + + Used for heat of vaporization and surface tension. + + Parameters + ---------- + Tc : float + Critical temperature :math:`a_0` [K]. + a1, a2, a3, a4, a5 : float + DIPPR-4 coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, Tc, a1, a2, a3=0.0, a4=0.0, a5=0.0): + self.coeffs = (Tc, a1, a2, a3, a4, a5) + super().__init__(func=self._eval) + + def _eval(self, T): + Tc, a1, a2, a3, a4, a5 = self.coeffs + Tr = T / Tc + h = a2 + a3 * Tr + a4 * Tr**2 + a5 * Tr**3 + return a1 * (1 - Tr)**h + + +class DIPPR5(Function): + r"""DIPPR Equation 5 correlation (IK-CAPE code: DIP5). + + .. math:: + + f(T) = \frac{a_0 \, T^{a_1}}{1 + a_2 / T + a_3 / T^2} + + Used for vapor viscosity and thermal conductivity. + + Parameters + ---------- + a0, a1, a2, a3 : float + DIPPR-5 coefficients. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"value": 0} + + def __init__(self, a0, a1, a2=0.0, a3=0.0): + self.coeffs = (a0, a1, a2, a3) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3 = self.coeffs + return a0 * T**a1 / (1 + a2 / T + a3 / T**2) diff --git a/tests/thermodynamics/__init__.py b/tests/thermodynamics/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/thermodynamics/test_correlations.py b/tests/thermodynamics/test_correlations.py new file mode 100644 index 0000000..e1e99af --- /dev/null +++ b/tests/thermodynamics/test_correlations.py @@ -0,0 +1,383 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.correlations.py' +## +## IK-CAPE Pure Component Property Correlations +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + Polynomial, + ExponentialPolynomial, + Watson, + Antoine, + ExtendedAntoine, + Kirchhoff, + ExtendedKirchhoff, + Sutherland, + Wagner, + LiquidHeatCapacity, + ExtendedLiquidHeatCapacity, + DynamicViscosity, + Rackett, + AlyLee, + DIPPR4, + DIPPR5, +) + + +# HELPERS ============================================================================== + +def eval_block(block, T): + """Set input T, update, and return output value.""" + block.inputs[0] = T + block.update(None) + return block.outputs[0] + + +# TESTS ================================================================================ + +class TestPolynomial(unittest.TestCase): + + def test_init(self): + P = Polynomial(coeffs=[1.0, 2.0, 3.0]) + self.assertEqual(P.coeffs, (1.0, 2.0, 3.0)) + + def test_constant(self): + P = Polynomial(coeffs=[42.0]) + self.assertAlmostEqual(eval_block(P, 300), 42.0) + + def test_linear(self): + P = Polynomial(coeffs=[10.0, 0.5]) + # f(300) = 10 + 0.5*300 = 160 + self.assertAlmostEqual(eval_block(P, 300), 160.0) + + def test_quadratic(self): + P = Polynomial(coeffs=[1.0, 0.0, 0.01]) + # f(10) = 1 + 0 + 0.01*100 = 2 + self.assertAlmostEqual(eval_block(P, 10), 2.0) + + +class TestExponentialPolynomial(unittest.TestCase): + + def test_init(self): + EP = ExponentialPolynomial(coeffs=[2.0, 0.001]) + self.assertEqual(EP.coeffs, (2.0, 0.001)) + + def test_constant_exponent(self): + EP = ExponentialPolynomial(coeffs=[3.0]) + # f(T) = 10^3 = 1000 + self.assertAlmostEqual(eval_block(EP, 300), 1000.0) + + def test_evaluation(self): + EP = ExponentialPolynomial(coeffs=[2.0, 0.001]) + T = 300 + expected = 10**(2.0 + 0.001 * 300) + self.assertAlmostEqual(eval_block(EP, T), expected, places=5) + + +class TestWatson(unittest.TestCase): + + def test_init(self): + W = Watson(a0=5.0e7, a1=0.38, a2=647.096) + self.assertEqual(W.coeffs, (5.0e7, 0.38, 647.096, 0.0)) + + def test_at_low_temp(self): + # Water-like: Hvap = a0 * (Tc - T)^a1 + W = Watson(a0=5.2053e7, a1=0.3199, a2=647.096) + result = eval_block(W, 373.15) + # Should be a positive value (heat of vaporization) + self.assertGreater(result, 0) + + def test_with_offset(self): + W = Watson(a0=1.0, a1=1.0, a2=500.0, a3=100.0) + # f(300) = 1.0*(500-300)^1 + 100 = 300 + self.assertAlmostEqual(eval_block(W, 300), 300.0) + + +class TestAntoine(unittest.TestCase): + + def test_init(self): + A = Antoine(a0=23.2256, a1=3835.18, a2=-45.343) + self.assertEqual(A.coeffs, (23.2256, 3835.18, -45.343)) + + def test_water_boiling_point(self): + # Water Antoine coefficients (ln, Pa, K) + # At 373.15 K, vapor pressure should be ~101325 Pa + A = Antoine(a0=23.2256, a1=3835.18, a2=-45.343) + P = eval_block(A, 373.15) + self.assertAlmostEqual(P, 101325, delta=2000) + + def test_water_at_room_temp(self): + A = Antoine(a0=23.2256, a1=3835.18, a2=-45.343) + P = eval_block(A, 298.15) + # Water vapor pressure at 25C ~ 3170 Pa + self.assertAlmostEqual(P, 3170, delta=200) + + +class TestExtendedAntoine(unittest.TestCase): + + def test_init_defaults(self): + EA = ExtendedAntoine(a0=10.0, a1=-3000.0, a2=-50.0) + self.assertEqual(EA.coeffs, (10.0, -3000.0, -50.0, 0.0, 0.0, 0.0, 0.0)) + + def test_reduces_to_antoine(self): + # When a3=a4=a5=0, should match Antoine (with sign convention a1 -> -a1) + A = Antoine(a0=23.2256, a1=3835.18, a2=-45.343) + EA = ExtendedAntoine(a0=23.2256, a1=-3835.18, a2=-45.343) + T = 373.15 + self.assertAlmostEqual(eval_block(A, T), eval_block(EA, T), places=5) + + +class TestKirchhoff(unittest.TestCase): + + def test_init(self): + K = Kirchhoff(a0=73.649, a1=7258.2, a2=7.3037) + self.assertEqual(K.coeffs, (73.649, 7258.2, 7.3037)) + + def test_evaluation(self): + # Verify numerical correctness: ln f(T) = a0 - a1/T - a2*ln(T) + K = Kirchhoff(a0=73.649, a1=7258.2, a2=7.3037) + T = 373.15 + expected = np.exp(73.649 - 7258.2 / T - 7.3037 * np.log(T)) + self.assertAlmostEqual(eval_block(K, T), expected, places=5) + + +class TestExtendedKirchhoff(unittest.TestCase): + + def test_init_defaults(self): + EK = ExtendedKirchhoff(a0=73.0, a1=-7000.0, a2=-7.0) + self.assertEqual(EK.coeffs, (73.0, -7000.0, -7.0, 0.0, 0.0)) + + def test_reduces_to_kirchhoff(self): + # When a3=0, should match Kirchhoff (with sign convention) + K = Kirchhoff(a0=73.649, a1=7258.2, a2=7.3037) + EK = ExtendedKirchhoff(a0=73.649, a1=-7258.2, a2=-7.3037) + T = 373.15 + self.assertAlmostEqual(eval_block(K, T), eval_block(EK, T), places=5) + + +class TestSutherland(unittest.TestCase): + + def test_init(self): + S = Sutherland(a0=1.458e-6, a1=110.4) + self.assertEqual(S.coeffs, (1.458e-6, 110.4)) + + def test_air_viscosity(self): + # Sutherland formula for air viscosity + # mu = a0*sqrt(T) / (1 + a1/T) + # At 300K, air viscosity ~ 1.85e-5 Pa.s + S = Sutherland(a0=1.458e-6, a1=110.4) + mu = eval_block(S, 300) + self.assertAlmostEqual(mu, 1.85e-5, delta=0.1e-5) + + def test_increases_with_temperature(self): + S = Sutherland(a0=1.458e-6, a1=110.4) + mu_300 = eval_block(S, 300) + mu_500 = eval_block(S, 500) + self.assertGreater(mu_500, mu_300) + + +class TestWagner(unittest.TestCase): + + def test_init(self): + W = Wagner(Tc=647.096, Pc=22064000, a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303) + self.assertEqual(W.coeffs, (647.096, 22064000, -7.76451, 1.45838, -2.77580, -1.23303)) + + def test_water_at_critical(self): + # At T = Tc, tau = 0, so f(Tc) = Pc + W = Wagner(Tc=647.096, Pc=22064000, a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303) + P = eval_block(W, 647.096) + self.assertAlmostEqual(P, 22064000, delta=1) + + def test_water_boiling_point(self): + # Wagner equation for water at 373.15K + W = Wagner(Tc=647.096, Pc=22064000, a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303) + P = eval_block(W, 373.15) + self.assertAlmostEqual(P, 101325, delta=2000) + + +class TestLiquidHeatCapacity(unittest.TestCase): + + def test_init_defaults(self): + C = LiquidHeatCapacity(a0=75.0) + self.assertEqual(C.coeffs, (75.0, 0.0, 0.0, 0.0, 0.0)) + + def test_constant_cp(self): + C = LiquidHeatCapacity(a0=75.4) + # f(T) = 75.4 when all other coeffs are zero + self.assertAlmostEqual(eval_block(C, 300), 75.4) + + def test_water_cp(self): + # Liquid water Cp ~ 75.3 J/(mol*K) at 298K, roughly constant + C = LiquidHeatCapacity(a0=75.3) + self.assertAlmostEqual(eval_block(C, 298.15), 75.3, places=1) + + +class TestExtendedLiquidHeatCapacity(unittest.TestCase): + + def test_init_defaults(self): + C = ExtendedLiquidHeatCapacity(a0=75.0) + self.assertEqual(C.coeffs, (75.0, 0.0, 0.0, 0.0, 0.0, 0.0)) + + def test_evaluation(self): + C = ExtendedLiquidHeatCapacity(a0=10.0, a1=0.1, a5=1000.0) + T = 200 + expected = 10.0 + 0.1 * 200 + 1000.0 / 200 + self.assertAlmostEqual(eval_block(C, T), expected) + + +class TestDynamicViscosity(unittest.TestCase): + + def test_init_default(self): + V = DynamicViscosity(a0=0.001, a1=1000) + self.assertEqual(V.coeffs, (0.001, 1000, 0.0)) + + def test_evaluation(self): + V = DynamicViscosity(a0=0.001, a1=1000) + T = 300 + expected = 0.001 * np.exp(1000 / 300) + self.assertAlmostEqual(eval_block(V, T), expected, places=8) + + def test_decreases_with_temperature(self): + # Liquid viscosity decreases with temperature + V = DynamicViscosity(a0=0.001, a1=1000) + mu_300 = eval_block(V, 300) + mu_400 = eval_block(V, 400) + self.assertGreater(mu_300, mu_400) + + +class TestRackett(unittest.TestCase): + + def test_init(self): + R = Rackett(a0=0.01805, a1=0.2882, a2=647.096, a3=0.2857) + self.assertEqual(R.coeffs, (0.01805, 0.2882, 647.096, 0.2857)) + + def test_density_positive(self): + R = Rackett(a0=0.01805, a1=0.2882, a2=647.096, a3=0.2857) + result = eval_block(R, 373.15) + self.assertGreater(result, 0) + + def test_density_decreases_with_temp(self): + # Liquid density decreases as T increases toward Tc + R = Rackett(a0=0.01805, a1=0.2882, a2=647.096, a3=0.2857) + rho_300 = eval_block(R, 300) + rho_500 = eval_block(R, 500) + self.assertGreater(rho_300, rho_500) + + +class TestAlyLee(unittest.TestCase): + + def test_init(self): + AL = AlyLee(a0=33363, a1=26790, a2=2610.5, a3=8896, a4=1169) + self.assertEqual(AL.coeffs, (33363, 26790, 2610.5, 8896, 1169)) + + def test_water_ideal_gas_cp(self): + # DIPPR Aly-Lee coefficients for water ideal gas Cp (J/kmol/K) + AL = AlyLee(a0=33363, a1=26790, a2=2610.5, a3=8896, a4=1169) + Cp = eval_block(AL, 300) + # Water ideal gas Cp at 300K ~ 33580 J/kmol/K + self.assertAlmostEqual(Cp, 33580, delta=500) + + def test_increases_with_temperature(self): + AL = AlyLee(a0=33363, a1=26790, a2=2610.5, a3=8896, a4=1169) + Cp_300 = eval_block(AL, 300) + Cp_500 = eval_block(AL, 500) + self.assertGreater(Cp_500, Cp_300) + + +class TestDIPPR4(unittest.TestCase): + + def test_init_defaults(self): + D = DIPPR4(Tc=647.096, a1=5.2053e7, a2=0.3199) + self.assertEqual(D.coeffs, (647.096, 5.2053e7, 0.3199, 0.0, 0.0, 0.0)) + + def test_zero_at_tc(self): + # At T = Tc, (1-Tr)=0, so f(Tc)=0 + D = DIPPR4(Tc=647.096, a1=5.2053e7, a2=0.3199) + self.assertAlmostEqual(eval_block(D, 647.096), 0.0) + + def test_water_hvap(self): + # Water heat of vaporization via DIPPR-4 + # At 373.15K, Hvap ~ 2.257e6 J/kg = ~40.65 kJ/mol + D = DIPPR4(Tc=647.096, a1=5.2053e7, a2=0.3199) + Hvap = eval_block(D, 373.15) + # Should be roughly 40.65e6 mJ/kmol range - ~4e7 J/kmol + self.assertGreater(Hvap, 0) + self.assertAlmostEqual(Hvap, 4.07e7, delta=0.5e7) + + +class TestDIPPR5(unittest.TestCase): + + def test_init_defaults(self): + D = DIPPR5(a0=1.0e-6, a1=0.5) + self.assertEqual(D.coeffs, (1.0e-6, 0.5, 0.0, 0.0)) + + def test_evaluation(self): + D = DIPPR5(a0=2.0, a1=1.5, a2=100.0, a3=5000.0) + T = 400 + expected = 2.0 * 400**1.5 / (1 + 100 / 400 + 5000 / 400**2) + self.assertAlmostEqual(eval_block(D, T), expected, places=5) + + def test_simple_power_law(self): + # When a2=a3=0, reduces to a0*T^a1 + D = DIPPR5(a0=3.0, a1=2.0) + T = 10 + self.assertAlmostEqual(eval_block(D, T), 300.0) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllBlocks(unittest.TestCase): + """Verify every block can be instantiated and evaluated at T=300K.""" + + def test_all_blocks_at_300K(self): + blocks = [ + Polynomial(coeffs=[1.0, 0.01]), + ExponentialPolynomial(coeffs=[2.0, 0.001]), + Watson(a0=5e7, a1=0.38, a2=647.0), + Antoine(a0=23.2256, a1=3835.18, a2=-45.343), + ExtendedAntoine(a0=23.2256, a1=-3835.18, a2=-45.343), + Kirchhoff(a0=73.649, a1=7258.2, a2=7.3037), + ExtendedKirchhoff(a0=73.649, a1=-7258.2, a2=-7.3037), + Sutherland(a0=1.458e-6, a1=110.4), + Wagner(Tc=647.096, Pc=22064000, a2=-7.76, a3=1.46, a4=-2.78, a5=-1.23), + LiquidHeatCapacity(a0=75.3), + ExtendedLiquidHeatCapacity(a0=75.3), + DynamicViscosity(a0=0.001, a1=1000), + Rackett(a0=0.01805, a1=0.2882, a2=647.096, a3=0.2857), + AlyLee(a0=33363, a1=26790, a2=2610.5, a3=8896, a4=1169), + DIPPR4(Tc=647.096, a1=5.2053e7, a2=0.3199), + DIPPR5(a0=1e-6, a1=0.5), + ] + + for block in blocks: + with self.subTest(block=block.__class__.__name__): + result = eval_block(block, 300) + self.assertTrue(np.isfinite(result), f"{block.__class__.__name__} returned non-finite: {result}") + + +# PORT LABELS ========================================================================== + +class TestPortLabels(unittest.TestCase): + """Verify all blocks have correct port labels.""" + + def test_input_port_label(self): + block = Antoine(a0=23.0, a1=3800.0, a2=-45.0) + self.assertEqual(block.input_port_labels, {"T": 0}) + + def test_output_port_label(self): + block = Antoine(a0=23.0, a1=3800.0, a2=-45.0) + self.assertEqual(block.output_port_labels, {"value": 0}) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From 2cd7e66363d321448cbf1660996d3a69d7e587de Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 11:08:58 +0100 Subject: [PATCH 2/9] Add IK-CAPE calculation of averages (10 mixing rule blocks) --- src/pathsim_chem/thermodynamics/__init__.py | 1 + src/pathsim_chem/thermodynamics/averages.py | 330 ++++++++++++++++++++ tests/thermodynamics/test_averages.py | 256 +++++++++++++++ 3 files changed, 587 insertions(+) create mode 100644 src/pathsim_chem/thermodynamics/averages.py create mode 100644 tests/thermodynamics/test_averages.py diff --git a/src/pathsim_chem/thermodynamics/__init__.py b/src/pathsim_chem/thermodynamics/__init__.py index f4a6307..7bcda50 100644 --- a/src/pathsim_chem/thermodynamics/__init__.py +++ b/src/pathsim_chem/thermodynamics/__init__.py @@ -1 +1,2 @@ from .correlations import * +from .averages import * diff --git a/src/pathsim_chem/thermodynamics/averages.py b/src/pathsim_chem/thermodynamics/averages.py new file mode 100644 index 0000000..2ee9c27 --- /dev/null +++ b/src/pathsim_chem/thermodynamics/averages.py @@ -0,0 +1,330 @@ +######################################################################################### +## +## IK-CAPE Calculation of Averages (Chapter 3) +## +## Mixing rules for computing mixture properties from pure component values +## and composition data. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# BLOCKS ================================================================================ + +class MolarAverage(Function): + r"""Molar average mixing rule (IK-CAPE code: MOLA). + + .. math:: + + \text{average} = \sum_i x_i \, \text{value}_i + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x): + self.x = np.asarray(x, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + return np.dot(self.x, values) + + +class MassAverage(Function): + r"""Mass fraction average mixing rule (IK-CAPE code: MASS). + + .. math:: + + \text{average} = \frac{\sum_i w_i \, \text{value}_i}{\sum_i w_i} + + Parameters + ---------- + w : array_like + Mass fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, w): + self.w = np.asarray(w, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + return np.dot(self.w, values) / np.sum(self.w) + + +class LogMolarAverage(Function): + r"""Logarithmic molar average mixing rule (IK-CAPE code: MOLG). + + .. math:: + + \ln(\text{average}) = \sum_i x_i \, \ln(\text{value}_i) + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x): + self.x = np.asarray(x, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + return np.exp(np.dot(self.x, np.log(values))) + + +class LogMassAverage(Function): + r"""Logarithmic mass fraction average mixing rule (IK-CAPE code: MALG). + + .. math:: + + \ln(\text{average}) = \frac{\sum_i w_i \, \ln(\text{value}_i)}{\sum_i w_i} + + Parameters + ---------- + w : array_like + Mass fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, w): + self.w = np.asarray(w, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + return np.exp(np.dot(self.w, np.log(values)) / np.sum(self.w)) + + +class LambdaAverage(Function): + r"""Thermal conductivity average for gaseous mixtures (IK-CAPE code: LAMB). + + .. math:: + + \lambda_m = 0.5 \left( \sum_i x_i \, \lambda_i + + \frac{1}{\sum_i x_i / \lambda_i} \right) + + Arithmetic-harmonic mean of pure component conductivities. + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x): + self.x = np.asarray(x, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + v = np.asarray(values) + return 0.5 * (np.dot(self.x, v) + 1.0 / np.dot(self.x, 1.0 / v)) + + +class ViscosityAverage(Function): + r"""Viscosity average for gaseous mixtures (IK-CAPE code: VISC). + + .. math:: + + \mu_m = \frac{\sum_i x_i \, \sqrt{M_i} \, \mu_i} + {\sum_i x_i \, \sqrt{M_i}} + + Molecular-weight-weighted average. + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + M : array_like + Molecular weights [kg/kmol] for each component. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x, M): + self.x = np.asarray(x, dtype=float) + self.M = np.asarray(M, dtype=float) + self.weights = self.x * np.sqrt(self.M) + super().__init__(func=self._eval) + + def _eval(self, *values): + return np.dot(self.weights, values) / np.sum(self.weights) + + +class VolumeAverage(Function): + r"""Volume-based density average (IK-CAPE code: VOLU). + + .. math:: + + \rho_m = \frac{1}{\sum_i x_i / \rho_i} + + Harmonic mean weighted by mole fractions, used for mixture density. + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x): + self.x = np.asarray(x, dtype=float) + super().__init__(func=self._eval) + + def _eval(self, *values): + return 1.0 / np.dot(self.x, 1.0 / np.asarray(values)) + + +class WilkeViscosity(Function): + r"""Wilke mixing rule for gas viscosity (IK-CAPE code: WILK). + + .. math:: + + \mu_m = \sum_i \frac{y_i \, \mu_i}{\sum_j y_j \, F_{ij}} + + .. math:: + + F_{ij} = \frac{\left(1 + \sqrt{\mu_i / \mu_j} + \, \sqrt[4]{M_j / M_i}\right)^2} + {\sqrt{8 \left(1 + M_i / M_j\right)}} + + Parameters + ---------- + y : array_like + Mole fractions (vapor phase) for each component. + M : array_like + Molecular weights [kg/kmol] for each component. + """ + + output_port_labels = {"average": 0} + + def __init__(self, y, M): + self.y = np.asarray(y, dtype=float) + self.M = np.asarray(M, dtype=float) + n = len(self.M) + + # precompute the molecular weight ratios + self._Mr = np.zeros((n, n)) + self._denom = np.zeros((n, n)) + for i in range(n): + for j in range(n): + self._Mr[i, j] = (self.M[j] / self.M[i])**0.25 + self._denom[i, j] = np.sqrt(8 * (1 + self.M[i] / self.M[j])) + + super().__init__(func=self._eval) + + def _eval(self, *values): + v = np.asarray(values) + n = len(v) + result = 0.0 + for i in range(n): + denom_sum = 0.0 + for j in range(n): + F_ij = (1 + np.sqrt(v[i] / v[j]) * self._Mr[i, j])**2 / self._denom[i, j] + denom_sum += self.y[j] * F_ij + result += self.y[i] * v[i] / denom_sum + return result + + +class WassiljewaMasonSaxena(Function): + r"""Wassiljewa-Mason-Saxena mixing rule for gas thermal conductivity + (IK-CAPE code: WAMA). + + .. math:: + + \lambda_m = \sum_i \frac{y_i \, \lambda_i}{\sum_j y_j \, F_{ij}} + + where :math:`F_{ij}` uses viscosity values :math:`\eta`: + + .. math:: + + F_{ij} = \frac{\left(1 + \sqrt{\eta_i / \eta_j} + \, \sqrt[4]{M_j / M_i}\right)^2} + {\sqrt{8 \left(1 + M_i / M_j\right)}} + + The first N inputs are thermal conductivities, the next N are viscosities. + + Parameters + ---------- + y : array_like + Mole fractions (vapor phase) for each component. + M : array_like + Molecular weights [kg/kmol] for each component. + """ + + output_port_labels = {"average": 0} + + def __init__(self, y, M): + self.y = np.asarray(y, dtype=float) + self.M = np.asarray(M, dtype=float) + self.n = len(self.M) + + # precompute molecular weight ratios + n = self.n + self._Mr = np.zeros((n, n)) + self._denom = np.zeros((n, n)) + for i in range(n): + for j in range(n): + self._Mr[i, j] = (self.M[j] / self.M[i])**0.25 + self._denom[i, j] = np.sqrt(8 * (1 + self.M[i] / self.M[j])) + + super().__init__(func=self._eval) + + def _eval(self, *inputs): + n = self.n + lam = np.asarray(inputs[:n]) # thermal conductivities + eta = np.asarray(inputs[n:]) # viscosities + + result = 0.0 + for i in range(n): + denom_sum = 0.0 + for j in range(n): + F_ij = (1 + np.sqrt(eta[i] / eta[j]) * self._Mr[i, j])**2 / self._denom[i, j] + denom_sum += self.y[j] * F_ij + result += self.y[i] * lam[i] / denom_sum + return result + + +class DIPPRSurfaceTension(Function): + r"""DIPPR surface tension average (IK-CAPE code: DIST). + + .. math:: + + \text{average} = \left( + \frac{\sum_i x_i \, V_i \, \text{value}_i^{1/4}} + {\sum_i x_i \, V_i} \right)^4 + + Parameters + ---------- + x : array_like + Mole fractions for each component. Determines the number of input ports. + V : array_like + Molar volumes [m^3/kmol] for each component. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x, V): + self.x = np.asarray(x, dtype=float) + self.V = np.asarray(V, dtype=float) + self.xV = self.x * self.V + super().__init__(func=self._eval) + + def _eval(self, *values): + v = np.asarray(values) + return (np.dot(self.xV, v**0.25) / np.sum(self.xV))**4 diff --git a/tests/thermodynamics/test_averages.py b/tests/thermodynamics/test_averages.py new file mode 100644 index 0000000..4e641c6 --- /dev/null +++ b/tests/thermodynamics/test_averages.py @@ -0,0 +1,256 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.averages.py' +## +## IK-CAPE Calculation of Averages (Chapter 3) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + MolarAverage, + MassAverage, + LogMolarAverage, + LogMassAverage, + LambdaAverage, + ViscosityAverage, + VolumeAverage, + WilkeViscosity, + WassiljewaMasonSaxena, + DIPPRSurfaceTension, +) + + +# HELPERS ============================================================================== + +def eval_block(block, *values): + """Set inputs and return output.""" + for i, v in enumerate(values): + block.inputs[i] = v + block.update(None) + return block.outputs[0] + + +# TESTS ================================================================================ + +class TestMolarAverage(unittest.TestCase): + + def test_init(self): + M = MolarAverage(x=[0.5, 0.5]) + np.testing.assert_array_equal(M.x, [0.5, 0.5]) + + def test_equal_fractions(self): + M = MolarAverage(x=[0.5, 0.5]) + # avg = 0.5*10 + 0.5*20 = 15 + self.assertAlmostEqual(eval_block(M, 10, 20), 15.0) + + def test_pure_component(self): + M = MolarAverage(x=[1.0, 0.0]) + self.assertAlmostEqual(eval_block(M, 42, 99), 42.0) + + def test_three_components(self): + M = MolarAverage(x=[0.2, 0.3, 0.5]) + # 0.2*100 + 0.3*200 + 0.5*300 = 20 + 60 + 150 = 230 + self.assertAlmostEqual(eval_block(M, 100, 200, 300), 230.0) + + +class TestMassAverage(unittest.TestCase): + + def test_normalized_fractions(self): + M = MassAverage(w=[0.4, 0.6]) + # avg = (0.4*10 + 0.6*20) / (0.4+0.6) = 16 + self.assertAlmostEqual(eval_block(M, 10, 20), 16.0) + + def test_unnormalized_fractions(self): + # w=[2, 3] should give same result as [0.4, 0.6] + M = MassAverage(w=[2, 3]) + self.assertAlmostEqual(eval_block(M, 10, 20), 16.0) + + +class TestLogMolarAverage(unittest.TestCase): + + def test_equal_values(self): + L = LogMolarAverage(x=[0.5, 0.5]) + # ln(avg) = 0.5*ln(100) + 0.5*ln(100) => avg = 100 + self.assertAlmostEqual(eval_block(L, 100, 100), 100.0) + + def test_geometric_mean(self): + L = LogMolarAverage(x=[0.5, 0.5]) + # geometric mean of 4 and 9 = sqrt(36) = 6 + self.assertAlmostEqual(eval_block(L, 4, 9), 6.0) + + def test_pure_component(self): + L = LogMolarAverage(x=[1.0, 0.0]) + self.assertAlmostEqual(eval_block(L, 42, 99), 42.0) + + +class TestLogMassAverage(unittest.TestCase): + + def test_geometric_mean(self): + L = LogMassAverage(w=[0.5, 0.5]) + # same as LogMolar when w sums to 1 + self.assertAlmostEqual(eval_block(L, 4, 9), 6.0) + + def test_unnormalized(self): + L = LogMassAverage(w=[1, 1]) + self.assertAlmostEqual(eval_block(L, 4, 9), 6.0) + + +class TestLambdaAverage(unittest.TestCase): + + def test_equal_values(self): + L = LambdaAverage(x=[0.5, 0.5]) + # 0.5*(sum + 1/sum_inv) with equal values = value + self.assertAlmostEqual(eval_block(L, 10, 10), 10.0) + + def test_arithmetic_harmonic_mean(self): + L = LambdaAverage(x=[0.5, 0.5]) + v1, v2 = 2.0, 8.0 + arith = 0.5 * (0.5 * v1 + 0.5 * v2) + harmo = 0.5 * (1.0 / (0.5 / v1 + 0.5 / v2)) + expected = arith + harmo + self.assertAlmostEqual(eval_block(L, v1, v2), expected) + + def test_pure_component(self): + L = LambdaAverage(x=[1.0, 0.0, 0.0]) + self.assertAlmostEqual(eval_block(L, 5.0, 10.0, 20.0), 5.0) + + +class TestViscosityAverage(unittest.TestCase): + + def test_equal_mol_weights(self): + # With equal M, reduces to molar average + V = ViscosityAverage(x=[0.5, 0.5], M=[28, 28]) + self.assertAlmostEqual(eval_block(V, 10, 20), 15.0) + + def test_weighting(self): + V = ViscosityAverage(x=[0.5, 0.5], M=[4, 16]) + # weights: 0.5*sqrt(4)=1, 0.5*sqrt(16)=2, sum=3 + # avg = (1*10 + 2*20) / 3 = 50/3 + self.assertAlmostEqual(eval_block(V, 10, 20), 50.0 / 3.0) + + +class TestVolumeAverage(unittest.TestCase): + + def test_equal_values(self): + V = VolumeAverage(x=[0.5, 0.5]) + self.assertAlmostEqual(eval_block(V, 1000, 1000), 1000.0) + + def test_harmonic_mean(self): + V = VolumeAverage(x=[0.5, 0.5]) + # 1 / (0.5/800 + 0.5/1200) = 1 / (0.000625 + 0.000417) = 960 + expected = 1.0 / (0.5 / 800 + 0.5 / 1200) + self.assertAlmostEqual(eval_block(V, 800, 1200), expected, places=5) + + def test_pure_component(self): + V = VolumeAverage(x=[1.0, 0.0]) + self.assertAlmostEqual(eval_block(V, 800, 1200), 800.0) + + +class TestWilkeViscosity(unittest.TestCase): + + def test_pure_component(self): + W = WilkeViscosity(y=[1.0, 0.0], M=[28, 32]) + # Pure component 1, should return its viscosity + self.assertAlmostEqual(eval_block(W, 1.8e-5, 2.0e-5), 1.8e-5) + + def test_equal_species(self): + # Two identical species => average should equal their shared value + W = WilkeViscosity(y=[0.5, 0.5], M=[28, 28]) + self.assertAlmostEqual(eval_block(W, 1.5e-5, 1.5e-5), 1.5e-5) + + def test_symmetric_F(self): + # For equal mol weights, F_ij depends only on viscosity ratio + W = WilkeViscosity(y=[0.5, 0.5], M=[28, 28]) + result = eval_block(W, 1.0e-5, 2.0e-5) + # Should be between the two values + self.assertGreater(result, 1.0e-5) + self.assertLess(result, 2.0e-5) + + def test_n2_o2_air(self): + # N2/O2 mixture (roughly air) + W = WilkeViscosity(y=[0.79, 0.21], M=[28.014, 31.999]) + mu_N2 = 1.78e-5 # Pa.s at ~300K + mu_O2 = 2.06e-5 + result = eval_block(W, mu_N2, mu_O2) + # Air viscosity ~ 1.85e-5 Pa.s + self.assertAlmostEqual(result, 1.85e-5, delta=0.1e-5) + + +class TestWassiljewaMasonSaxena(unittest.TestCase): + + def test_pure_component(self): + W = WassiljewaMasonSaxena(y=[1.0, 0.0], M=[28, 32]) + # inputs: lambda_1, lambda_2, eta_1, eta_2 + lam1, lam2 = 0.026, 0.027 + eta1, eta2 = 1.8e-5, 2.0e-5 + result = eval_block(W, lam1, lam2, eta1, eta2) + self.assertAlmostEqual(result, lam1) + + def test_equal_species(self): + W = WassiljewaMasonSaxena(y=[0.5, 0.5], M=[28, 28]) + lam, eta = 0.025, 1.5e-5 + result = eval_block(W, lam, lam, eta, eta) + self.assertAlmostEqual(result, lam) + + def test_bounded(self): + W = WassiljewaMasonSaxena(y=[0.5, 0.5], M=[28, 32]) + lam1, lam2 = 0.020, 0.030 + eta1, eta2 = 1.5e-5, 2.0e-5 + result = eval_block(W, lam1, lam2, eta1, eta2) + self.assertGreater(result, lam1) + self.assertLess(result, lam2) + + +class TestDIPPRSurfaceTension(unittest.TestCase): + + def test_equal_values(self): + D = DIPPRSurfaceTension(x=[0.5, 0.5], V=[0.05, 0.05]) + # When all values are equal, average = value + self.assertAlmostEqual(eval_block(D, 0.072, 0.072), 0.072, places=6) + + def test_pure_component(self): + D = DIPPRSurfaceTension(x=[1.0, 0.0], V=[0.05, 0.04]) + self.assertAlmostEqual(eval_block(D, 0.072, 0.030), 0.072, places=6) + + def test_bounded(self): + D = DIPPRSurfaceTension(x=[0.5, 0.5], V=[0.05, 0.05]) + result = eval_block(D, 0.030, 0.072) + self.assertGreater(result, 0.030) + self.assertLess(result, 0.072) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllAverages(unittest.TestCase): + """Verify every averaging block can be instantiated and evaluated.""" + + def test_all_blocks(self): + blocks_and_inputs = [ + (MolarAverage(x=[0.5, 0.5]), [10, 20]), + (MassAverage(w=[0.4, 0.6]), [10, 20]), + (LogMolarAverage(x=[0.5, 0.5]), [10, 20]), + (LogMassAverage(w=[0.5, 0.5]), [10, 20]), + (LambdaAverage(x=[0.5, 0.5]), [10, 20]), + (ViscosityAverage(x=[0.5, 0.5], M=[28, 32]), [1e-5, 2e-5]), + (VolumeAverage(x=[0.5, 0.5]), [800, 1200]), + (WilkeViscosity(y=[0.5, 0.5], M=[28, 32]), [1.5e-5, 2e-5]), + (WassiljewaMasonSaxena(y=[0.5, 0.5], M=[28, 32]), [0.025, 0.026, 1.5e-5, 2e-5]), + (DIPPRSurfaceTension(x=[0.5, 0.5], V=[0.05, 0.04]), [0.03, 0.07]), + ] + + for block, inputs in blocks_and_inputs: + with self.subTest(block=block.__class__.__name__): + result = eval_block(block, *inputs) + self.assertTrue(np.isfinite(result), f"{block.__class__.__name__} returned non-finite: {result}") + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From 64e3e9b29df52fc1ab0ce3813754f7021d8c20fb Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 14:31:44 +0100 Subject: [PATCH 3/9] Add activity coefficient models (NRTL, Wilson, UNIQUAC) and DIKL average --- src/pathsim_chem/thermodynamics/__init__.py | 3 + .../thermodynamics/activity_coefficients.py | 334 ++++++++++++++++++ src/pathsim_chem/thermodynamics/averages.py | 37 ++ .../test_activity_coefficients.py | 242 +++++++++++++ 4 files changed, 616 insertions(+) create mode 100644 src/pathsim_chem/thermodynamics/activity_coefficients.py create mode 100644 tests/thermodynamics/test_activity_coefficients.py diff --git a/src/pathsim_chem/thermodynamics/__init__.py b/src/pathsim_chem/thermodynamics/__init__.py index 7bcda50..5649713 100644 --- a/src/pathsim_chem/thermodynamics/__init__.py +++ b/src/pathsim_chem/thermodynamics/__init__.py @@ -1,2 +1,5 @@ from .correlations import * from .averages import * +from .activity_coefficients import * +from .equations_of_state import * +from .corrections import * diff --git a/src/pathsim_chem/thermodynamics/activity_coefficients.py b/src/pathsim_chem/thermodynamics/activity_coefficients.py new file mode 100644 index 0000000..90cd4a0 --- /dev/null +++ b/src/pathsim_chem/thermodynamics/activity_coefficients.py @@ -0,0 +1,334 @@ +######################################################################################### +## +## IK-CAPE Activity Coefficient Models (Chapter 4) +## +## Models for computing liquid-phase activity coefficients from +## composition and temperature. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# BLOCKS ================================================================================ + +class NRTL(Function): + r"""Non-Random Two-Liquid activity coefficient model (IK-CAPE Chapter 4.1). + + Computes liquid-phase activity coefficients for a multi-component mixture + at a given temperature. The NRTL model is widely used for strongly + non-ideal liquid mixtures, including partially miscible systems. + It is especially suitable for polar and associating systems such as + alcohol-water, amine-water, and organic acid mixtures. + + The model uses binary interaction parameters :math:`\tau_{ij}` and + non-randomness parameters :math:`\alpha_{ij}` that are typically fitted + to experimental VLE data. + + **Input port:** ``T`` -- temperature [K]. + + **Output ports:** ``out_0``, ``out_1``, ... -- activity coefficients + :math:`\gamma_i` for each component (one per mole fraction entry). + + .. math:: + + \ln \gamma_i = \frac{\sum_j \tau_{ji} G_{ji} x_j}{\sum_k G_{ki} x_k} + + \sum_j \frac{x_j G_{ij}}{\sum_k G_{kj} x_k} + \left( \tau_{ij} - \frac{\sum_l \tau_{lj} G_{lj} x_l}{\sum_k G_{kj} x_k} \right) + + where: + + .. math:: + + G_{ji} = \exp(-\alpha_{ji} \, \tau_{ji}) + + \tau_{ji} = a_{ji} + b_{ji}/T + e_{ji} \ln T + f_{ji} T + + \alpha_{ji} = c_{ji} + d_{ji} (T - 273.15) + + Parameters + ---------- + x : array_like + Mole fractions [N]. Fixed mixture composition. + a : array_like + Constant part of the energy interaction parameter :math:`\tau` [N x N]. + Diagonal elements should be zero. For constant (temperature-independent) + tau values, set ``a=tau`` and leave ``b``, ``e``, ``f`` as zero. + b : array_like, optional + Coefficient of :math:`1/T` in :math:`\tau_{ji}` [N x N]. Default: zeros. + c : array_like, optional + Constant part of the non-randomness parameter :math:`\alpha` [N x N]. + Default: 0.3 for all off-diagonal pairs (a common starting value). + d : array_like, optional + Temperature-dependent part of :math:`\alpha` [N x N]. Default: zeros. + e : array_like, optional + Coefficient of :math:`\ln T` in :math:`\tau_{ji}` [N x N]. Default: zeros. + f : array_like, optional + Coefficient of :math:`T` in :math:`\tau_{ji}` [N x N]. Default: zeros. + """ + + input_port_labels = {"T": 0} + + def __init__(self, x, a, b=None, c=None, d=None, e=None, f=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.e = np.zeros((n, n)) if e is None else np.asarray(e, dtype=float).reshape(n, n) + self.f = np.zeros((n, n)) if f is None else np.asarray(f, dtype=float).reshape(n, n) + + if c is None: + self.c = np.full((n, n), 0.3) + np.fill_diagonal(self.c, 0.0) + else: + self.c = np.asarray(c, dtype=float).reshape(n, n) + + self.d = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # temperature-dependent parameters + tau = self.a + self.b / T + self.e * np.log(T) + self.f * T + alpha = self.c + self.d * (T - 273.15) + G = np.exp(-alpha * tau) + + ln_gamma = np.zeros(n) + for i in range(n): + # denominators: sum_k G_ki * x_k for each column + den_i = np.dot(G[:, i], x) + + # first term + num1 = np.dot(tau[:, i] * G[:, i], x) + term1 = num1 / den_i + + # second term + term2 = 0.0 + for j in range(n): + den_j = np.dot(G[:, j], x) + num_inner = np.dot(tau[:, j] * G[:, j], x) + term2 += (x[j] * G[i, j] / den_j) * (tau[i, j] - num_inner / den_j) + + ln_gamma[i] = term1 + term2 + + return tuple(np.exp(ln_gamma)) + + +class Wilson(Function): + r"""Wilson activity coefficient model (IK-CAPE Chapter 4.3). + + Computes liquid-phase activity coefficients for a multi-component mixture. + The Wilson model is well suited for completely miscible systems with + moderate non-ideality, such as hydrocarbon mixtures and many + organic-organic systems. It cannot represent liquid-liquid phase splitting. + + The model uses binary interaction parameters :math:`\Lambda_{ij}` related + to molar volume ratios and energy differences between molecular pairs. + + **Input port:** ``T`` -- temperature [K]. + + **Output ports:** ``out_0``, ``out_1``, ... -- activity coefficients + :math:`\gamma_i` for each component. + + .. math:: + + \ln \gamma_i = 1 - \ln\!\left(\sum_j x_j \Lambda_{ij}\right) + - \sum_k \frac{x_k \Lambda_{ki}}{\sum_j x_j \Lambda_{kj}} + + where: + + .. math:: + + \Lambda_{ij} = \exp(a_{ij} + b_{ij}/T + c_{ij} \ln T + d_{ij} T) + + Parameters + ---------- + x : array_like + Mole fractions [N]. Fixed mixture composition. + a : array_like + Constant part of :math:`\ln \Lambda_{ij}` [N x N]. Diagonal should + be zero so that :math:`\Lambda_{ii} = 1`. + b : array_like, optional + Coefficient of :math:`1/T` in :math:`\ln \Lambda_{ij}` [N x N]. + Default: zeros. + c : array_like, optional + Coefficient of :math:`\ln T` in :math:`\ln \Lambda_{ij}` [N x N]. + Default: zeros. + d : array_like, optional + Coefficient of :math:`T` in :math:`\ln \Lambda_{ij}` [N x N]. + Default: zeros. + """ + + input_port_labels = {"T": 0} + + def __init__(self, x, a, b=None, c=None, d=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.c = np.zeros((n, n)) if c is None else np.asarray(c, dtype=float).reshape(n, n) + self.d = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # temperature-dependent Lambda + Lam = np.exp(self.a + self.b / T + self.c * np.log(T) + self.d * T) + + ln_gamma = np.zeros(n) + for i in range(n): + # sum_j x_j * Lambda_ij + s1 = np.dot(x, Lam[i, :]) + + # sum_k x_k * Lambda_ki / sum_j x_j * Lambda_kj + s2 = 0.0 + for k in range(n): + den_k = np.dot(x, Lam[k, :]) + s2 += x[k] * Lam[k, i] / den_k + + ln_gamma[i] = 1.0 - np.log(s1) - s2 + + return tuple(np.exp(ln_gamma)) + + +class UNIQUAC(Function): + r"""Universal Quasi-Chemical activity coefficient model (IK-CAPE Chapter 4.2). + + Computes liquid-phase activity coefficients by combining a combinatorial + contribution (based on molecular size and shape via van der Waals volume + ``r`` and surface area ``q`` parameters) with a residual contribution + (based on intermolecular interactions via temperature-dependent + :math:`\tau_{ij}` parameters). UNIQUAC can handle strongly non-ideal + systems including partially miscible liquids and forms the theoretical + basis of the UNIFAC group-contribution method. + + **Input port:** ``T`` -- temperature [K]. + + **Output ports:** ``out_0``, ``out_1``, ... -- activity coefficients + :math:`\gamma_i` for each component. + + The activity coefficient is split into combinatorial and residual parts: + + .. math:: + + \ln \gamma_i = \ln \gamma_i^C + \ln \gamma_i^R + + Combinatorial part (molecular size/shape effects): + + .. math:: + + \ln \gamma_i^C = \ln\frac{V_i}{x_i} + \frac{z}{2} q_i \ln\frac{F_i}{V_i} + + l_i - \frac{V_i}{x_i} \sum_j x_j l_j + + Residual part (intermolecular energy interactions): + + .. math:: + + \ln \gamma_i^R = q'_i \left[1 - \ln\!\left(\sum_j F'_j \tau_{ji}\right) + - \sum_j \frac{F'_j \tau_{ij}}{\sum_k F'_k \tau_{kj}} \right] + + Parameters + ---------- + x : array_like + Mole fractions [N]. Fixed mixture composition. + r : array_like + Van der Waals volume (size) parameters [N], from UNIFAC tables + or fitted to data. Example: ethanol=2.1055, water=0.92. + q : array_like + Van der Waals surface area parameters [N]. Example: ethanol=1.972, + water=1.4. + a : array_like + Constant part of the interaction parameter exponent [N x N]. + Diagonal should be zero. The interaction parameter is computed as + :math:`\tau_{ji} = \exp(a_{ji} + b_{ji}/T + c_{ji} \ln T + d_{ji} T)`. + q_prime : array_like, optional + Modified surface area for the residual part [N]. Used in some + formulations (e.g., for water/alcohol systems). Defaults to ``q``. + b : array_like, optional + Coefficient of :math:`1/T` in the :math:`\tau` exponent [N x N]. + Default: zeros. + c : array_like, optional + Coefficient of :math:`\ln T` in the :math:`\tau` exponent [N x N]. + Default: zeros. + d : array_like, optional + Coefficient of :math:`T` in the :math:`\tau` exponent [N x N]. + Default: zeros. + z : float, optional + Lattice coordination number (default 10, the standard value). + """ + + input_port_labels = {"T": 0} + + def __init__(self, x, r, q, a, q_prime=None, b=None, c=None, d=None, z=10): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.r = np.asarray(r, dtype=float) + self.q = np.asarray(q, dtype=float) + self.q_prime = self.q.copy() if q_prime is None else np.asarray(q_prime, dtype=float) + self.z = z + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b_param = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.c_param = np.zeros((n, n)) if c is None else np.asarray(c, dtype=float).reshape(n, n) + self.d_param = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + # l_i = z/2 * (r_i - q_i) - (r_i - 1) + self.l = self.z / 2 * (self.r - self.q) - (self.r - 1) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # temperature-dependent interaction parameters + 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) + Fp = self.q_prime * x / np.dot(self.q_prime, x) + + # combinatorial part + 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) + + # residual part + ln_gamma_R = np.zeros(n) + for i in range(n): + s1 = np.dot(Fp, tau[:, i]) # sum_j F'_j * tau_ji + + s2 = 0.0 + for j in range(n): + den = np.dot(Fp, tau[:, j]) # sum_k F'_k * tau_kj + s2 += Fp[j] * tau[i, j] / den + + ln_gamma_R[i] = self.q_prime[i] * (1 - np.log(s1) - s2) + + ln_gamma = ln_gamma_C + ln_gamma_R + return tuple(np.exp(ln_gamma)) diff --git a/src/pathsim_chem/thermodynamics/averages.py b/src/pathsim_chem/thermodynamics/averages.py index 2ee9c27..ea82c0a 100644 --- a/src/pathsim_chem/thermodynamics/averages.py +++ b/src/pathsim_chem/thermodynamics/averages.py @@ -328,3 +328,40 @@ def __init__(self, x, V): def _eval(self, *values): v = np.asarray(values) return (np.dot(self.xV, v**0.25) / np.sum(self.xV))**4 + + +class DIPPRLiquidConductivity(Function): + r"""DIPPR liquid thermal conductivity average (IK-CAPE code: DIKL). + + .. math:: + + \text{average} = \sum_i \sum_j \frac{2 \, xv_i \, xv_j} + {1/\text{value}_i + 1/\text{value}_j} + + where :math:`xv_i = x_i V_i / \sum_j x_j V_j` are the volume fractions. + + Parameters + ---------- + x : array_like + Mole fractions for each component. + V : array_like + Molar volumes [m^3/kmol] for each component. + """ + + output_port_labels = {"average": 0} + + def __init__(self, x, V): + self.x = np.asarray(x, dtype=float) + self.V = np.asarray(V, dtype=float) + xV = self.x * self.V + self.xv = xV / np.sum(xV) + super().__init__(func=self._eval) + + def _eval(self, *values): + v = np.asarray(values) + n = len(v) + result = 0.0 + for i in range(n): + for j in range(n): + result += 2 * self.xv[i] * self.xv[j] / (1 / v[i] + 1 / v[j]) + return result diff --git a/tests/thermodynamics/test_activity_coefficients.py b/tests/thermodynamics/test_activity_coefficients.py new file mode 100644 index 0000000..43c687c --- /dev/null +++ b/tests/thermodynamics/test_activity_coefficients.py @@ -0,0 +1,242 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.activity_coefficients.py' +## +## IK-CAPE Activity Coefficient Models (Chapter 4) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + NRTL, + Wilson, + UNIQUAC, +) + + +# HELPERS ============================================================================== + +def eval_block(block, T): + """Set input T, update, and return all outputs.""" + block.inputs[0] = T + block.update(None) + n_out = len([k for k in block.outputs]) + if n_out == 1: + return block.outputs[0] + return tuple(block.outputs[i] for i in range(n_out)) + + +# TESTS ================================================================================ + +class TestNRTL(unittest.TestCase): + + def test_init(self): + N = NRTL( + x=[0.5, 0.5], + a=[[0, 1.5], [2.0, 0]], + ) + self.assertEqual(N.n, 2) + np.testing.assert_array_equal(N.x, [0.5, 0.5]) + + def test_pure_component_gamma_is_one(self): + # For x_1=1.0, gamma_1 should be 1.0 regardless of parameters + N = NRTL( + x=[1.0, 0.0], + a=[[0, 1.5], [2.0, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + gamma1, gamma2 = eval_block(N, 350) + self.assertAlmostEqual(gamma1, 1.0, places=10) + + def test_symmetric_mixture(self): + # Symmetric parameters and equal x => gamma_1 = gamma_2 + N = NRTL( + x=[0.5, 0.5], + a=[[0, 1.0], [1.0, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + gamma1, gamma2 = eval_block(N, 350) + self.assertAlmostEqual(gamma1, gamma2, places=10) + self.assertGreater(gamma1, 1.0) # positive deviation + + def test_ethanol_water(self): + # Ethanol(1)-Water(2) NRTL parameters (simplified, constant tau and alpha) + # tau_12 = -0.801, tau_21 = 3.458, alpha = 0.3 + N = NRTL( + x=[0.4, 0.6], + a=[[0, -0.801], [3.458, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + gamma1, gamma2 = eval_block(N, 350) + # Both should be > 1 for this system + self.assertGreater(gamma1, 1.0) + self.assertGreater(gamma2, 1.0) + + def test_temperature_dependence(self): + # With b parameter, gamma should change with T + N = NRTL( + x=[0.5, 0.5], + a=[[0, 0.5], [0.5, 0]], + b=[[0, 100], [100, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + g300 = eval_block(N, 300) + g400 = eval_block(N, 400) + # gamma should change with temperature + self.assertNotAlmostEqual(g300[0], g400[0], places=3) + + def test_three_components(self): + N = NRTL( + x=[0.3, 0.3, 0.4], + a=[[0, 0.5, 0.8], + [0.6, 0, 0.4], + [0.9, 0.3, 0]], + c=[[0, 0.3, 0.3], + [0.3, 0, 0.3], + [0.3, 0.3, 0]], + ) + gammas = eval_block(N, 350) + self.assertEqual(len(gammas), 3) + for g in gammas: + self.assertGreater(g, 0) + self.assertTrue(np.isfinite(g)) + + +class TestWilson(unittest.TestCase): + + def test_init(self): + W = Wilson( + x=[0.5, 0.5], + a=[[0, 0.5], [-0.3, 0]], + ) + self.assertEqual(W.n, 2) + + def test_pure_component_gamma_is_one(self): + W = Wilson( + x=[1.0, 0.0], + a=[[0, 0.5], [-0.3, 0]], + ) + gamma1, gamma2 = eval_block(W, 350) + self.assertAlmostEqual(gamma1, 1.0, places=10) + + def test_symmetric_mixture(self): + W = Wilson( + x=[0.5, 0.5], + a=[[0, 0.5], [0.5, 0]], + ) + gamma1, gamma2 = eval_block(W, 350) + self.assertAlmostEqual(gamma1, gamma2, places=10) + + def test_identity_lambdas(self): + # When a=0 (all Lambdas = 1), ideal solution => gamma = 1 + W = Wilson( + x=[0.5, 0.5], + a=[[0, 0], [0, 0]], + ) + gamma1, gamma2 = eval_block(W, 350) + self.assertAlmostEqual(gamma1, 1.0, places=10) + self.assertAlmostEqual(gamma2, 1.0, places=10) + + def test_temperature_dependence(self): + W = Wilson( + x=[0.5, 0.5], + a=[[0, 0.2], [-0.3, 0]], + b=[[0, 200], [-150, 0]], + ) + g300 = eval_block(W, 300) + g400 = eval_block(W, 400) + self.assertNotAlmostEqual(g300[0], g400[0], places=3) + + +class TestUNIQUAC(unittest.TestCase): + + def test_init(self): + U = UNIQUAC( + x=[0.5, 0.5], + r=[2.1055, 0.92], + q=[1.972, 1.4], + a=[[0, -1.318], [2.772, 0]], + ) + self.assertEqual(U.n, 2) + + def test_pure_component_gamma_is_one(self): + U = UNIQUAC( + x=[1.0, 0.0], + r=[2.1055, 0.92], + q=[1.972, 1.4], + a=[[0, -1.318], [2.772, 0]], + ) + gamma1, gamma2 = eval_block(U, 350) + self.assertAlmostEqual(gamma1, 1.0, places=10) + + def test_symmetric_case(self): + # Identical r, q, symmetric a => gamma_1 = gamma_2 + U = UNIQUAC( + x=[0.5, 0.5], + r=[1.0, 1.0], + q=[1.0, 1.0], + a=[[0, -0.5], [-0.5, 0]], + ) + gamma1, gamma2 = eval_block(U, 350) + self.assertAlmostEqual(gamma1, gamma2, places=10) + + def test_ethanol_water(self): + # Ethanol(1)-Water(2) UNIQUAC + # r: ethanol=2.1055, water=0.92 + # q: ethanol=1.972, water=1.4 + U = UNIQUAC( + x=[0.4, 0.6], + r=[2.1055, 0.92], + q=[1.972, 1.4], + a=[[0, -1.318], [2.772, 0]], + ) + gamma1, gamma2 = eval_block(U, 350) + self.assertGreater(gamma1, 0) + self.assertGreater(gamma2, 0) + self.assertTrue(np.isfinite(gamma1)) + self.assertTrue(np.isfinite(gamma2)) + + def test_three_components(self): + U = UNIQUAC( + x=[0.3, 0.3, 0.4], + r=[2.1, 0.92, 1.5], + q=[1.97, 1.4, 1.2], + a=[[0, -1.0, 0.5], + [2.0, 0, -0.3], + [0.8, 0.4, 0]], + ) + gammas = eval_block(U, 350) + self.assertEqual(len(gammas), 3) + for g in gammas: + self.assertGreater(g, 0) + self.assertTrue(np.isfinite(g)) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllActivityCoefficients(unittest.TestCase): + + def test_all_at_350K(self): + blocks = [ + NRTL(x=[0.5, 0.5], a=[[0, 1.0], [1.5, 0]], c=[[0, 0.3], [0.3, 0]]), + Wilson(x=[0.5, 0.5], a=[[0, 0.5], [-0.3, 0]]), + UNIQUAC(x=[0.5, 0.5], r=[2.1, 0.92], q=[1.97, 1.4], a=[[0, -1.0], [2.0, 0]]), + ] + + for block in blocks: + with self.subTest(block=block.__class__.__name__): + result = eval_block(block, 350) + for g in result: + self.assertTrue(np.isfinite(g), f"{block.__class__.__name__} non-finite") + self.assertGreater(g, 0) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From 61d1657083b9f5dc4c41bef94fc796ff53596af7 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 14:31:52 +0100 Subject: [PATCH 4/9] Add equations of state (Peng-Robinson, SRK) and corrections (Poynting, Henry) --- .../thermodynamics/corrections.py | 115 +++++++ .../thermodynamics/equations_of_state.py | 309 ++++++++++++++++++ tests/thermodynamics/test_corrections.py | 115 +++++++ .../thermodynamics/test_equations_of_state.py | 169 ++++++++++ 4 files changed, 708 insertions(+) create mode 100644 src/pathsim_chem/thermodynamics/corrections.py create mode 100644 src/pathsim_chem/thermodynamics/equations_of_state.py create mode 100644 tests/thermodynamics/test_corrections.py create mode 100644 tests/thermodynamics/test_equations_of_state.py diff --git a/src/pathsim_chem/thermodynamics/corrections.py b/src/pathsim_chem/thermodynamics/corrections.py new file mode 100644 index 0000000..68d96cb --- /dev/null +++ b/src/pathsim_chem/thermodynamics/corrections.py @@ -0,0 +1,115 @@ +######################################################################################### +## +## IK-CAPE Corrections (Chapters 5 and 6) +## +## Poynting correction and Henry's law constant for vapor-liquid +## equilibrium calculations. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# BLOCKS ================================================================================ + +class PoyntingCorrection(Function): + r"""Poynting pressure correction factor (IK-CAPE Chapter 5). + + Accounts for the effect of total system pressure on the fugacity of a + liquid component when the system pressure :math:`P` differs from the + pure component saturation pressure :math:`P_i^s`. The Poynting factor + is used in rigorous VLE calculations via the modified Raoult's law: + + .. math:: + + y_i P = x_i \gamma_i P_i^s \phi_i^s F_{p,i} / \phi_i^V + + At moderate pressures (below ~10 bar), the Poynting correction is + very close to unity and can often be neglected. It becomes significant + at high pressures or for components with large liquid molar volumes. + + .. math:: + + \ln F_{p,i} = \frac{v_i^L \, (P - P_i^s)}{R T} + + **Input ports:** + + - ``T`` -- temperature [K] + - ``P`` -- total system pressure [Pa] + - ``Psat`` -- pure component saturation pressure [Pa] (e.g., from Antoine block) + - ``vL`` -- liquid molar volume [m^3/mol] (e.g., from Rackett block) + + **Output port:** ``Fp`` -- dimensionless Poynting correction factor [-]. + + Parameters + ---------- + None. All values are received as dynamic inputs from upstream blocks. + """ + + input_port_labels = {"T": 0, "P": 1, "Psat": 2, "vL": 3} + output_port_labels = {"Fp": 0} + + def __init__(self): + super().__init__(func=self._eval) + + def _eval(self, T, P, Psat, vL): + return np.exp(vL * (P - Psat) / (R * T)) + + +class HenryConstant(Function): + r"""Temperature-dependent Henry's law constant (IK-CAPE Chapter 6). + + Computes the Henry's law constant :math:`H_{i,j}` for a dissolved gas + species *i* in solvent *j* as a function of temperature. Henry's law + relates the partial pressure of a dilute gas above a solution to its + mole fraction in the liquid: + + .. math:: + + p_i = H_{i,j} \, x_i + + The temperature dependence follows the four-parameter correlation: + + .. math:: + + \ln H_{i,j} = a + \frac{b}{T} + c \, \ln T + d \, T + + Coefficients are specific to each gas-solvent pair and are available + from standard databases (e.g., DECHEMA, NIST). + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``H`` -- Henry's law constant [Pa]. + + Parameters + ---------- + a : float + Constant term in :math:`\ln H`. + b : float, optional + Coefficient of :math:`1/T`. Controls the temperature sensitivity + (related to enthalpy of dissolution). Default: 0. + c : float, optional + Coefficient of :math:`\ln T`. Default: 0. + d : float, optional + Coefficient of :math:`T`. Default: 0. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"H": 0} + + def __init__(self, a, b=0.0, c=0.0, d=0.0): + self.coeffs = (a, b, c, d) + super().__init__(func=self._eval) + + def _eval(self, T): + a, b, c, d = self.coeffs + return np.exp(a + b / T + c * np.log(T) + d * T) diff --git a/src/pathsim_chem/thermodynamics/equations_of_state.py b/src/pathsim_chem/thermodynamics/equations_of_state.py new file mode 100644 index 0000000..e4508d8 --- /dev/null +++ b/src/pathsim_chem/thermodynamics/equations_of_state.py @@ -0,0 +1,309 @@ +######################################################################################### +## +## IK-CAPE Equations of State (Chapter 7) +## +## Cubic equations of state for computing molar volumes and +## compressibility factors of gas and liquid phases. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# HELPERS =============================================================================== + +def _solve_cubic_eos(coeffs, phase="vapor"): + """Solve cubic equation in Z and return the appropriate root. + + Parameters + ---------- + coeffs : array_like + Polynomial coefficients [1, c2, c1, c0] for Z^3 + c2*Z^2 + c1*Z + c0 = 0. + phase : str + "vapor" selects the largest real root, "liquid" the smallest positive real root. + + Returns + ------- + float + The selected compressibility factor Z. + """ + roots = np.roots(coeffs) + + # keep only real roots (small imaginary part) + real_roots = roots[np.abs(roots.imag) < 1e-10].real + + # filter positive roots + positive_roots = real_roots[real_roots > 0] + + if len(positive_roots) == 0: + # fallback: take root with smallest imaginary part + idx = np.argmin(np.abs(roots.imag)) + return roots[idx].real + + if phase == "vapor": + return float(np.max(positive_roots)) + else: + return float(np.min(positive_roots)) + + +# BLOCKS ================================================================================ + +class PengRobinson(Function): + r"""Peng-Robinson cubic equation of state (IK-CAPE Chapter 7.2). + + Computes the molar volume and compressibility factor of a pure component + or mixture at a given temperature and pressure. The Peng-Robinson EoS + (1976) is one of the most widely used cubic equations of state in + chemical engineering, offering improved liquid density predictions + over the original Redlich-Kwong equation. It is suitable for + hydrocarbons, gases, and many industrial fluids. + + Supports both pure components (single Tc/Pc/omega) and mixtures + (arrays of Tc/Pc/omega plus mole fractions and optional binary + interaction parameters). The cubic equation is solved analytically + and the ``phase`` parameter controls whether the vapor (largest Z) + or liquid (smallest Z) root is selected. + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output ports:** ``v`` -- molar volume [m^3/mol], + ``z`` -- compressibility factor [-]. + + Pressure-explicit form: + + .. math:: + + P = \frac{RT}{v - b} - \frac{a(T)}{v(v+b) + b(v-b)} + + Solved as a cubic in compressibility factor :math:`Z = Pv/(RT)`: + + .. math:: + + Z^3 - (1-B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0 + + where :math:`A = a_m P / (R^2 T^2)` and :math:`B = b_m P / (RT)`. + + Pure component parameters are computed from critical properties: + + .. math:: + + a_i = 0.45724 \frac{R^2 T_{c,i}^2}{P_{c,i}} \alpha_i(T), \quad + b_i = 0.07780 \frac{R T_{c,i}}{P_{c,i}} + + \alpha_i = \left[1 + m_i(1 - \sqrt{T/T_{c,i}})\right]^2, \quad + m_i = 0.37464 + 1.54226\omega_i - 0.26992\omega_i^2 + + Standard van der Waals one-fluid mixing rules: + + .. math:: + + a_m = \sum_i \sum_j x_i x_j \sqrt{a_i a_j}(1 - k_{ij}), \quad + b_m = \sum_i x_i b_i + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. Scalar for pure component, array for mixture. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. Characterizes deviation from spherical + symmetry. Available in standard reference tables (e.g., DIPPR). + x : array_like, optional + Mole fractions [N]. Required for mixtures, omit for pure components. + k : array_like, optional + Binary interaction parameters [N x N]. Symmetric, with + :math:`k_{ii} = 0`. Default: zero for all pairs. + phase : str, optional + Phase root selection: ``"vapor"`` (default) picks the largest Z root, + ``"liquid"`` picks the smallest positive Z root. + """ + + input_port_labels = {"T": 0, "P": 1} + output_port_labels = {"v": 0, "z": 1} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.n = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.n == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.n, self.n)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.n, self.n) + + # constant pure component parameters + self.m = 0.37464 + 1.54226 * self.omega - 0.26992 * self.omega**2 + self.a_c = 0.45724 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.07780 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + n = self.n + x = self.x + + # temperature-dependent alpha and a_i + alpha = (1 + self.m * (1 - np.sqrt(T / self.Tc)))**2 + a_i = self.a_c * alpha + + # mixing rules + if n == 1: + a_m = a_i[0] + b_m = self.b_i[0] + else: + a_m = 0.0 + for i in range(n): + for j in range(n): + a_ij = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + a_m += x[i] * x[j] * a_ij + b_m = np.dot(x, self.b_i) + + # dimensionless parameters + A = a_m * P / (R**2 * T**2) + B = b_m * P / (R * T) + + # cubic in Z: Z^3 - (1-B)*Z^2 + (A-3B^2-2B)*Z - (AB-B^2-B^3) = 0 + coeffs = [1, + -(1 - B), + A - 3 * B**2 - 2 * B, + -(A * B - B**2 - B**3)] + + Z = _solve_cubic_eos(coeffs, self.phase) + v = Z * R * T / P + + return v, Z + + +class RedlichKwongSoave(Function): + r"""Soave-Redlich-Kwong cubic equation of state (IK-CAPE Chapter 7.1). + + Computes the molar volume and compressibility factor of a pure component + or mixture at a given temperature and pressure. The SRK EoS (Soave, 1972) + introduced a temperature-dependent attractive term to the original + Redlich-Kwong equation, making it suitable for VLE calculations of + hydrocarbons and light gases. It is widely used in natural gas + and refinery process simulation. + + Has the same interface and mixing rules as :class:`PengRobinson`, but + uses different constants and a slightly different pressure equation + (the attractive term denominator is :math:`v(v+b)` instead of + :math:`v(v+b)+b(v-b)`). + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output ports:** ``v`` -- molar volume [m^3/mol], + ``z`` -- compressibility factor [-]. + + Pressure-explicit form: + + .. math:: + + P = \frac{RT}{v - b} - \frac{a(T)}{v(v + b)} + + Solved as a cubic in :math:`Z`: + + .. math:: + + Z^3 - Z^2 + (A - B - B^2)Z - AB = 0 + + Pure component parameters: + + .. math:: + + a_i = 0.42748 \frac{R^2 T_{c,i}^2}{P_{c,i}} \alpha_i(T), \quad + b_i = 0.08664 \frac{R T_{c,i}}{P_{c,i}} + + \alpha_i = \left[1 + m_i(1 - \sqrt{T/T_{c,i}})\right]^2, \quad + m_i = 0.48 + 1.574\omega_i - 0.176\omega_i^2 + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. Scalar for pure component, array for mixture. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. + x : array_like, optional + Mole fractions [N]. Required for mixtures, omit for pure components. + k : array_like, optional + Binary interaction parameters [N x N]. Default: zero for all pairs. + phase : str, optional + Phase root selection: ``"vapor"`` (default) or ``"liquid"``. + """ + + input_port_labels = {"T": 0, "P": 1} + output_port_labels = {"v": 0, "z": 1} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.n = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.n == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.n, self.n)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.n, self.n) + + # constant pure component parameters + self.m = 0.48 + 1.574 * self.omega - 0.176 * self.omega**2 + self.a_c = 0.42748 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.08664 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + n = self.n + x = self.x + + # temperature-dependent alpha and a_i + alpha = (1 + self.m * (1 - np.sqrt(T / self.Tc)))**2 + a_i = self.a_c * alpha + + # mixing rules + if n == 1: + a_m = a_i[0] + b_m = self.b_i[0] + else: + a_m = 0.0 + for i in range(n): + for j in range(n): + a_ij = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + a_m += x[i] * x[j] * a_ij + b_m = np.dot(x, self.b_i) + + # dimensionless parameters + A = a_m * P / (R**2 * T**2) + B = b_m * P / (R * T) + + # cubic in Z: Z^3 - Z^2 + (A-B-B^2)*Z - AB = 0 + coeffs = [1, -1, A - B - B**2, -A * B] + + Z = _solve_cubic_eos(coeffs, self.phase) + v = Z * R * T / P + + return v, Z diff --git a/tests/thermodynamics/test_corrections.py b/tests/thermodynamics/test_corrections.py new file mode 100644 index 0000000..0b05d04 --- /dev/null +++ b/tests/thermodynamics/test_corrections.py @@ -0,0 +1,115 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.corrections.py' +## +## IK-CAPE Poynting Correction (Chapter 5) and Henry Constant (Chapter 6) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + PoyntingCorrection, + HenryConstant, +) + + +# CONSTANTS ============================================================================ + +R = 8.314462 + + +# HELPERS ============================================================================== + +def eval_block(block, *values): + """Set inputs and return output.""" + for i, v in enumerate(values): + block.inputs[i] = v + block.update(None) + return block.outputs[0] + + +# TESTS ================================================================================ + +class TestPoyntingCorrection(unittest.TestCase): + + def test_no_pressure_difference(self): + # When P == Psat, correction should be 1.0 + PC = PoyntingCorrection() + Fp = eval_block(PC, 373.15, 101325, 101325, 1.8e-5) + self.assertAlmostEqual(Fp, 1.0, places=10) + + def test_high_pressure(self): + # At P > Psat, Fp > 1 + PC = PoyntingCorrection() + Fp = eval_block(PC, 373.15, 1e7, 101325, 1.8e-5) + self.assertGreater(Fp, 1.0) + + def test_known_value(self): + # Manual calculation: Fp = exp(vL*(P-Psat)/(RT)) + T = 300 + P = 1e7 + Psat = 3500 + vL = 1.8e-5 # m^3/mol, typical liquid water + expected = np.exp(vL * (P - Psat) / (R * T)) + + PC = PoyntingCorrection() + Fp = eval_block(PC, T, P, Psat, vL) + self.assertAlmostEqual(Fp, expected, places=10) + + def test_close_to_one_at_moderate_pressure(self): + # For typical liquid, Poynting correction is small + PC = PoyntingCorrection() + Fp = eval_block(PC, 300, 200000, 100000, 1.8e-5) + # vL*(P-Psat)/(RT) ~ 1.8e-5 * 1e5 / (8.314*300) ~ 0.0007 + self.assertAlmostEqual(Fp, 1.0, delta=0.001) + + +class TestHenryConstant(unittest.TestCase): + + def test_init(self): + H = HenryConstant(a=15.0, b=-1500.0) + self.assertEqual(H.coeffs, (15.0, -1500.0, 0.0, 0.0)) + + def test_constant(self): + # With only a, H = exp(a) + H = HenryConstant(a=10.0) + result = eval_block(H, 300) + self.assertAlmostEqual(result, np.exp(10.0), places=5) + + def test_temperature_dependence(self): + H = HenryConstant(a=15.0, b=-1500.0) + H_300 = eval_block(H, 300) + H_350 = eval_block(H, 350) + # Different temperatures should give different values + self.assertNotAlmostEqual(H_300, H_350, places=0) + # Both should be positive + self.assertGreater(H_300, 0) + self.assertGreater(H_350, 0) + + def test_known_value(self): + # ln(H) = a + b/T + c*ln(T) + d*T + a, b, c, d = 10.0, -2000.0, 1.5, -0.001 + T = 350 + expected = np.exp(a + b / T + c * np.log(T) + d * T) + + H = HenryConstant(a=a, b=b, c=c, d=d) + result = eval_block(H, T) + self.assertAlmostEqual(result, expected, places=5) + + def test_o2_in_water(self): + # O2 Henry constant in water: ln(H/Pa) ~ 21 - 1700/T + H = HenryConstant(a=21.0, b=-1700.0) + H_val = eval_block(H, 298.15) + # Should be in the range of ~4e4 to ~5e4 kPa for O2 in water + self.assertGreater(H_val, 1e6) # positive and large + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) diff --git a/tests/thermodynamics/test_equations_of_state.py b/tests/thermodynamics/test_equations_of_state.py new file mode 100644 index 0000000..fe48c07 --- /dev/null +++ b/tests/thermodynamics/test_equations_of_state.py @@ -0,0 +1,169 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.equations_of_state.py' +## +## IK-CAPE Equations of State (Chapter 7) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + PengRobinson, + RedlichKwongSoave, +) + + +# CONSTANTS ============================================================================ + +R = 8.314462 + + +# HELPERS ============================================================================== + +def eval_block(block, T, P): + """Set T and P inputs, update, return (v, z).""" + block.inputs[0] = T + block.inputs[1] = P + block.update(None) + return block.outputs[0], block.outputs[1] + + +# TESTS ================================================================================ + +class TestPengRobinson(unittest.TestCase): + + def test_init_pure(self): + # Methane: Tc=190.6K, Pc=4.6MPa, omega=0.011 + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(PR.n, 1) + + def test_ideal_gas_limit(self): + # At low pressure and high T, Z -> 1 and v -> RT/P + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011) + T, P = 1000, 100 # very high T, very low P + v, z = eval_block(PR, T, P) + self.assertAlmostEqual(z, 1.0, delta=0.01) + expected_v = R * T / P + self.assertAlmostEqual(v, expected_v, delta=expected_v * 0.01) + + def test_methane_gas(self): + # Methane at 300K, 1 atm - should behave as near-ideal gas + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011) + v, z = eval_block(PR, 300, 101325) + self.assertAlmostEqual(z, 1.0, delta=0.05) + self.assertGreater(v, 0) + + def test_water_vapor(self): + # Water at 500K, 1 atm + PR = PengRobinson(Tc=647.096, Pc=22064000, omega=0.3443) + v, z = eval_block(PR, 500, 101325) + self.assertAlmostEqual(z, 1.0, delta=0.05) + + def test_liquid_phase(self): + # Methane below Tc at high pressure - liquid phase + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid") + v, z = eval_block(PR, 120, 5e6) + self.assertLess(z, 0.5) + self.assertGreater(v, 0) + + def test_mixture(self): + # Methane + ethane mixture + PR = PengRobinson( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + v, z = eval_block(PR, 300, 101325) + self.assertGreater(z, 0.9) + self.assertGreater(v, 0) + + def test_at_critical_point(self): + # At Tc, Pc => Z should be near 0.307 for PR + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011) + v, z = eval_block(PR, 190.6, 4.6e6) + self.assertAlmostEqual(z, 0.307, delta=0.05) + + +class TestRedlichKwongSoave(unittest.TestCase): + + def test_init_pure(self): + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(RKS.n, 1) + + def test_ideal_gas_limit(self): + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011) + T, P = 1000, 100 + v, z = eval_block(RKS, T, P) + self.assertAlmostEqual(z, 1.0, delta=0.01) + + def test_methane_gas(self): + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011) + v, z = eval_block(RKS, 300, 101325) + self.assertAlmostEqual(z, 1.0, delta=0.05) + self.assertGreater(v, 0) + + def test_liquid_phase(self): + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid") + v, z = eval_block(RKS, 120, 5e6) + self.assertLess(z, 0.5) + self.assertGreater(v, 0) + + def test_mixture(self): + RKS = RedlichKwongSoave( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + v, z = eval_block(RKS, 300, 101325) + self.assertGreater(z, 0.9) + self.assertGreater(v, 0) + + def test_pr_vs_rks_similar(self): + # Both EoS should give similar results for ideal gas conditions + PR = PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011) + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011) + + v_pr, z_pr = eval_block(PR, 300, 101325) + v_rks, z_rks = eval_block(RKS, 300, 101325) + + self.assertAlmostEqual(z_pr, z_rks, delta=0.01) + + def test_at_critical_point(self): + # At Tc, Pc => Z should be near 1/3 for RKS + RKS = RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011) + v, z = eval_block(RKS, 190.6, 4.6e6) + self.assertAlmostEqual(z, 0.333, delta=0.05) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllEoS(unittest.TestCase): + + def test_all_blocks(self): + blocks = [ + PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011), + PengRobinson(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid"), + RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011), + RedlichKwongSoave(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid"), + ] + + for block in blocks: + with self.subTest(block=f"{block.__class__.__name__}({block.phase})"): + v, z = eval_block(block, 200, 1e6) + self.assertTrue(np.isfinite(v)) + self.assertTrue(np.isfinite(z)) + self.assertGreater(v, 0) + self.assertGreater(z, 0) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From d991d981b58d8e3271cec311b7d100abdcb0c754 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 14:53:33 +0100 Subject: [PATCH 5/9] Add fugacity coefficients, enthalpy, reactions, and Flory-Huggins (Ch 4,8,9,10) --- src/pathsim_chem/thermodynamics/__init__.py | 3 + .../thermodynamics/activity_coefficients.py | 107 ++++ src/pathsim_chem/thermodynamics/enthalpy.py | 605 ++++++++++++++++++ .../thermodynamics/fugacity_coefficients.py | 353 ++++++++++ src/pathsim_chem/thermodynamics/reactions.py | 228 +++++++ .../test_activity_coefficients.py | 72 +++ tests/thermodynamics/test_enthalpy.py | 362 +++++++++++ .../test_fugacity_coefficients.py | 195 ++++++ tests/thermodynamics/test_reactions.py | 224 +++++++ 9 files changed, 2149 insertions(+) create mode 100644 src/pathsim_chem/thermodynamics/enthalpy.py create mode 100644 src/pathsim_chem/thermodynamics/fugacity_coefficients.py create mode 100644 src/pathsim_chem/thermodynamics/reactions.py create mode 100644 tests/thermodynamics/test_enthalpy.py create mode 100644 tests/thermodynamics/test_fugacity_coefficients.py create mode 100644 tests/thermodynamics/test_reactions.py diff --git a/src/pathsim_chem/thermodynamics/__init__.py b/src/pathsim_chem/thermodynamics/__init__.py index 5649713..b2ad04f 100644 --- a/src/pathsim_chem/thermodynamics/__init__.py +++ b/src/pathsim_chem/thermodynamics/__init__.py @@ -3,3 +3,6 @@ from .activity_coefficients import * from .equations_of_state import * from .corrections import * +from .fugacity_coefficients import * +from .enthalpy import * +from .reactions import * diff --git a/src/pathsim_chem/thermodynamics/activity_coefficients.py b/src/pathsim_chem/thermodynamics/activity_coefficients.py index 90cd4a0..a4ca198 100644 --- a/src/pathsim_chem/thermodynamics/activity_coefficients.py +++ b/src/pathsim_chem/thermodynamics/activity_coefficients.py @@ -332,3 +332,110 @@ def _eval(self, T): ln_gamma = ln_gamma_C + ln_gamma_R return tuple(np.exp(ln_gamma)) + + +class FloryHuggins(Function): + r"""Flory-Huggins activity coefficient model (IK-CAPE Chapter 4.4). + + Computes liquid-phase activity coefficients for polymer-solvent systems. + The Flory-Huggins model accounts for the large size difference between + polymer and solvent molecules by using volume fractions instead of mole + fractions. It is the standard model for polymer solutions and is suitable + for systems with one solvent and one or more dissolved polymer species. + + Component 1 is always the **solvent** (with segment number :math:`r_1 = 1`). + The remaining components are polymer species with segment numbers + :math:`r_{2i}` that represent the ratio of polymer to solvent molar volume. + + **Input port:** ``T`` -- temperature [K]. + + **Output ports:** ``out_0``, ``out_1``, ... -- activity coefficients + :math:`\gamma_i` for each component (solvent first, then polymers). + + The model uses volume fractions and the Flory-Huggins interaction + parameter :math:`\chi`: + + .. math:: + + \phi_1 = \frac{x_1}{\bar{r}}, \quad + \phi_{2i} = \frac{r_{2i} x_{2i}}{\bar{r}}, \quad + \bar{r} = x_1 + \sum_i r_{2i} x_{2i} + + Solvent activity coefficient: + + .. math:: + + \ln \gamma_1 = 1 + \ln\frac{1}{\bar{r}} - \frac{1}{\bar{r}} + + (1 - \phi_1) \sum_i \phi_{2i} \chi_{1,2i} + + Polymer activity coefficient for species k: + + .. math:: + + \ln \gamma_{2k} = 1 + \ln\frac{r_{2k}}{\bar{r}} - \frac{r_{2k}}{\bar{r}} + + \phi_1 r_{2k} \left( \chi_{1,2k} + - \sum_i \phi_{2i} \chi_{1,2i} \right) + + Temperature-dependent interaction parameter: + + .. math:: + + \chi_{1,2i} = \chi^0_{1,2i} + \frac{\chi^1_{1,2i}}{T} + + Parameters + ---------- + x : array_like + Mole fractions [N] where ``x[0]`` is the solvent and ``x[1:]`` + are the polymer species. + r : array_like + Segment numbers [N-1] for each polymer species. These represent + the ratio of the polymer molar volume to the solvent molar volume. + chi0 : array_like + Constant part of the Flory-Huggins interaction parameter [N-1]. + chi1 : array_like, optional + Temperature-dependent part of :math:`\chi` [N-1]. Default: zeros. + """ + + input_port_labels = {"T": 0} + + def __init__(self, x, r, chi0, chi1=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n_poly = self.n - 1 + + self.r_poly = np.asarray(r, dtype=float) + self.chi0 = np.asarray(chi0, dtype=float) + self.chi1 = np.zeros(n_poly) if chi1 is None else np.asarray(chi1, dtype=float) + + super().__init__(func=self._eval) + + def _eval(self, T): + x = self.x + x1 = x[0] + x2 = x[1:] + + # interaction parameters + chi = self.chi0 + self.chi1 / T + + # average segment number (r_solvent = 1) + r_bar = x1 + np.dot(self.r_poly, x2) + + # volume fractions + phi1 = x1 / r_bar + phi2 = self.r_poly * x2 / r_bar + + # weighted interaction sum + chi_sum = np.dot(phi2, chi) + + # solvent activity coefficient + ln_gamma1 = 1.0 + np.log(1.0 / r_bar) - 1.0 / r_bar + (1.0 - phi1) * chi_sum + + # polymer activity coefficients + ln_gamma2 = np.zeros(len(x2)) + for k in range(len(x2)): + ln_gamma2[k] = (1.0 + np.log(self.r_poly[k] / r_bar) + - self.r_poly[k] / r_bar + + phi1 * self.r_poly[k] * (chi[k] - chi_sum)) + + ln_gamma = np.concatenate(([ln_gamma1], ln_gamma2)) + return tuple(np.exp(ln_gamma)) diff --git a/src/pathsim_chem/thermodynamics/enthalpy.py b/src/pathsim_chem/thermodynamics/enthalpy.py new file mode 100644 index 0000000..24ff03e --- /dev/null +++ b/src/pathsim_chem/thermodynamics/enthalpy.py @@ -0,0 +1,605 @@ +######################################################################################### +## +## IK-CAPE Enthalpy (Chapter 9) +## +## Excess enthalpy models (9.6) and isothermal pressure dependency +## of the enthalpy in the vapor phase (9.7). +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + +from .equations_of_state import _solve_cubic_eos + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# EXCESS ENTHALPY MODELS (9.6) ========================================================== + +class ExcessEnthalpyNRTL(Function): + r"""Excess enthalpy from the NRTL model (IK-CAPE Chapter 9.6.1). + + Computes the molar excess enthalpy :math:`h^E` of a liquid mixture + using the NRTL (Non-Random Two-Liquid) model. The excess enthalpy + is obtained from the temperature derivative of the excess Gibbs energy: + + .. math:: + + h^E = -T^2 \frac{\partial (g^E / T)}{\partial T} + = -R T^2 \sum_i x_i \frac{A'_i B_i - A_i B'_i}{B_i^2} + + where: + + .. math:: + + A_i = \sum_j x_j G_{j,i} \tau_{j,i}, \quad + A'_i = \sum_j x_j (G'_{j,i} \tau_{j,i} + G_{j,i} \tau'_{j,i}) + + B_i = \sum_j x_j G_{j,i}, \quad + B'_i = \sum_j x_j G'_{j,i} + + Temperature derivatives: + + .. math:: + + \tau'_{j,i} = -b_{j,i}/T^2 + e_{j,i}/T + f_{j,i} + + G'_{j,i} = -G_{j,i}(d_{j,i} \tau_{j,i} + S_{j,i} \tau'_{j,i}) + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``hE`` -- molar excess enthalpy [J/mol]. + + Parameters + ---------- + x : array_like + Mole fractions [N]. Fixed mixture composition. + a : array_like + Constant part of :math:`\tau` [N x N]. Diagonal should be zero. + b : array_like, optional + Coefficient of :math:`1/T` in :math:`\tau` [N x N]. Default: zeros. + c : array_like, optional + Constant part of non-randomness :math:`\alpha` [N x N]. + Default: 0.3 off-diagonal. + d : array_like, optional + Temperature-dependent part of :math:`\alpha` [N x N]. Default: zeros. + e : array_like, optional + Coefficient of :math:`\ln T` in :math:`\tau` [N x N]. Default: zeros. + f : array_like, optional + Coefficient of :math:`T` in :math:`\tau` [N x N]. Default: zeros. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"hE": 0} + + def __init__(self, x, a, b=None, c=None, d=None, e=None, f=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.e = np.zeros((n, n)) if e is None else np.asarray(e, dtype=float).reshape(n, n) + self.f = np.zeros((n, n)) if f is None else np.asarray(f, dtype=float).reshape(n, n) + + if c is None: + self.c = np.full((n, n), 0.3) + np.fill_diagonal(self.c, 0.0) + else: + self.c = np.asarray(c, dtype=float).reshape(n, n) + + self.d = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # tau and its derivative + tau = self.a + self.b / T + self.e * np.log(T) + self.f * T + tau_prime = -self.b / T**2 + self.e / T + self.f + + # alpha (S) and its derivative + S = self.c + self.d * (T - 273.15) + S_prime = self.d + + # G and its derivative + G = np.exp(-S * tau) + G_prime = -G * (S_prime * tau + S * tau_prime) + + # h^E = -RT^2 * sum_i x_i * (A'_i * B_i - A_i * B'_i) / B_i^2 + hE = 0.0 + for i in range(n): + A_i = np.dot(x, G[:, i] * tau[:, i]) + A_prime_i = np.dot(x, G_prime[:, i] * tau[:, i] + G[:, i] * tau_prime[:, i]) + B_i = np.dot(x, G[:, i]) + B_prime_i = np.dot(x, G_prime[:, i]) + + hE += x[i] * (A_prime_i * B_i - A_i * B_prime_i) / B_i**2 + + hE = -R * T**2 * hE + return hE + + +class ExcessEnthalpyUNIQUAC(Function): + r"""Excess enthalpy from the UNIQUAC model (IK-CAPE Chapter 9.6.2). + + Computes the molar excess enthalpy :math:`h^E` of a liquid mixture + using the UNIQUAC model. Only the residual (interaction) part + contributes to the excess enthalpy since the combinatorial part + depends only on composition, not temperature: + + .. math:: + + h^E = R T^2 \sum_i q'_i x_i \left( + \frac{A'_i}{A_i} + \sum_j F'_j + \frac{\frac{\partial \tau_{i,j}}{\partial T} A_j + - \tau_{i,j} A'_j}{A_j^2} + \right) + + where :math:`A_i = \sum_j F'_j \tau_{j,i}` and + :math:`\tau_{j,i} = \exp(a_{j,i} + b_{j,i}/T + c_{j,i}\ln T + d_{j,i}T)`. + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``hE`` -- molar excess enthalpy [J/mol]. + + Parameters + ---------- + x : array_like + Mole fractions [N]. + r : array_like + Van der Waals volume parameters [N]. + q : array_like + Van der Waals surface area parameters [N]. + a : array_like + Constant part of the :math:`\\tau` exponent [N x N]. + q_prime : array_like, optional + Modified surface area for residual part [N]. Defaults to ``q``. + b : array_like, optional + Coefficient of :math:`1/T` [N x N]. Default: zeros. + c : array_like, optional + Coefficient of :math:`\\ln T` [N x N]. Default: zeros. + d : array_like, optional + Coefficient of :math:`T` [N x N]. Default: zeros. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"hE": 0} + + def __init__(self, x, r, q, a, q_prime=None, b=None, c=None, d=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.r = np.asarray(r, dtype=float) + self.q = np.asarray(q, dtype=float) + self.q_prime = self.q.copy() if q_prime is None else np.asarray(q_prime, dtype=float) + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b_param = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.c_param = np.zeros((n, n)) if c is None else np.asarray(c, dtype=float).reshape(n, n) + self.d_param = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # tau and derivative + exponent = self.a + self.b_param / T + self.c_param * np.log(T) + self.d_param * T + tau = np.exp(exponent) + tau_prime = tau * (-self.b_param / T**2 + self.c_param / T + self.d_param) + + # modified surface fractions F' + Fp = self.q_prime * x / np.dot(self.q_prime, x) + + # A_i = sum_j F'_j * tau_j,i and A'_i = sum_j F'_j * tau'_j,i + A_vec = np.dot(Fp, tau) # shape (n,) + A_prime = np.dot(Fp, tau_prime) # shape (n,) + + # h^E = RT^2 * sum_i q'_i * x_i * (A'_i/A_i + sum_j F'_j * (tau'_ij*A_j - tau_ij*A'_j)/A_j^2) + hE = 0.0 + for i in range(n): + term1 = A_prime[i] / A_vec[i] + + term2 = 0.0 + for j in range(n): + term2 += Fp[j] * (tau_prime[i, j] * A_vec[j] + - tau[i, j] * A_prime[j]) / A_vec[j]**2 + + hE += self.q_prime[i] * x[i] * (term1 + term2) + + return R * T**2 * hE + + +class ExcessEnthalpyWilson(Function): + r"""Excess enthalpy from the Wilson model (IK-CAPE Chapter 9.6.3). + + Computes the molar excess enthalpy :math:`h^E` of a liquid mixture + using the Wilson activity coefficient model. The Wilson model is + suitable for completely miscible systems: + + .. math:: + + \frac{g^E}{T} = -R \sum_i x_i \ln\!\left(\sum_j x_j \Lambda_{ij}\right) + + .. math:: + + h^E = R T^2 \sum_i x_i \frac{\sigma'_i}{\sigma_i} + + where :math:`\sigma_i = \sum_j x_j \Lambda_{ij}` and + :math:`\sigma'_i = \sum_j x_j \Lambda'_{ij}`, with: + + .. math:: + + \Lambda'_{ij} = \Lambda_{ij} + \left(-\frac{b_{ij}}{T^2} + \frac{c_{ij}}{T} + d_{ij}\right) + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``hE`` -- molar excess enthalpy [J/mol]. + + Parameters + ---------- + x : array_like + Mole fractions [N]. + a : array_like + Constant part of :math:`\ln \Lambda_{ij}` [N x N]. + b : array_like, optional + Coefficient of :math:`1/T` [N x N]. Default: zeros. + c : array_like, optional + Coefficient of :math:`\ln T` [N x N]. Default: zeros. + d : array_like, optional + Coefficient of :math:`T` [N x N]. Default: zeros. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"hE": 0} + + def __init__(self, x, a, b=None, c=None, d=None): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + n = self.n + + self.a = np.asarray(a, dtype=float).reshape(n, n) + self.b = np.zeros((n, n)) if b is None else np.asarray(b, dtype=float).reshape(n, n) + self.c = np.zeros((n, n)) if c is None else np.asarray(c, dtype=float).reshape(n, n) + self.d = np.zeros((n, n)) if d is None else np.asarray(d, dtype=float).reshape(n, n) + + super().__init__(func=self._eval) + + def _eval(self, T): + n = self.n + x = self.x + + # Lambda and its derivative + Lam = np.exp(self.a + self.b / T + self.c * np.log(T) + self.d * T) + Lam_prime = Lam * (-self.b / T**2 + self.c / T + self.d) + + # sigma_i = sum_j x_j * Lambda_ij, sigma'_i = sum_j x_j * Lambda'_ij + sigma = np.dot(Lam, x) # shape (n,) -- Lam[i,:] . x + sigma_prime = np.dot(Lam_prime, x) + + # h^E = RT^2 * sum_i x_i * sigma'_i / sigma_i + hE = R * T**2 * np.dot(x, sigma_prime / sigma) + return hE + + +class ExcessEnthalpyRedlichKister(Function): + r"""Excess enthalpy from the Redlich-Kister expansion (IK-CAPE Chapter 9.6.5). + + 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: + + .. math:: + + h^E = \frac{1}{2} \sum_i \sum_j h^E_{ij} + + Each binary pair contribution: + + .. 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) + + 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` + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``hE`` -- molar excess enthalpy [J/mol]. + + Parameters + ---------- + x : array_like + 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:: + + {(0, 1): [[a0, a1], [b0, b1], [c0]]} + + This gives A(T)=a0+a1*T, B(T)=b0+b1*T, C(T)=c0. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"hE": 0} + + def __init__(self, x, coeffs): + self.x = np.asarray(x, dtype=float) + self.n = len(self.x) + self.coeffs = coeffs + + super().__init__(func=self._eval) + + def _eval(self, T): + x = self.x + hE = 0.0 + + for (i, j), terms in self.coeffs.items(): + xi, xj = x[i], x[j] + x_sum = xi + xj + if x_sum < 1e-30: + continue + + xd = xi - xj + + # evaluate each Redlich-Kister term + pair_hE = 0.0 + xd_power = xd + for poly_coeffs in terms: + # temperature polynomial: c0 + c1*T + c2*T^2 + ... + coeff_val = 0.0 + T_power = 1.0 + for c in poly_coeffs: + coeff_val += c * T_power + T_power *= T + pair_hE += coeff_val * xd_power + xd_power *= xd + + hE += xi * xj / x_sum * pair_hE + + return 0.5 * hE + + +# ENTHALPY DEPARTURE FROM EOS (9.7) ===================================================== + +class EnthalpyDepartureRKS(Function): + r"""Isothermal enthalpy departure from the SRK EoS (IK-CAPE Chapter 9.7.1). + + Computes the departure of the real-gas enthalpy from the ideal-gas + enthalpy at the same temperature using the Soave-Redlich-Kwong equation + of state. This quantity is needed to convert between ideal-gas and + real-gas thermodynamic properties: + + .. math:: + + \Delta h = h(T,P) - h^{ig}(T) + + For a pure component: + + .. math:: + + \Delta h_i = RT(Z_i - 1) + \frac{T\frac{da_i}{dT} - a_i}{b_i} + \ln\!\left(1 + \frac{b_i P_i^s}{Z_i R T}\right) + + For a mixture: + + .. math:: + + \Delta h_m = RT(Z_m - 1) + \frac{T\frac{da_m}{dT} - a_m}{b_m} + \ln\!\left(1 + \frac{b_m P}{Z_m R T}\right) + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output port:** ``dh`` -- enthalpy departure [J/mol]. + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. + x : array_like, optional + Mole fractions [N]. Required for mixtures. + k : array_like, optional + Binary interaction parameters [N x N]. Default: zeros. + phase : str, optional + ``"vapor"`` (default) or ``"liquid"``. + """ + + input_port_labels = {"T": 0, "P": 1} + output_port_labels = {"dh": 0} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.nc = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.nc == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.nc, self.nc)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.nc, self.nc) + + # SRK constants + self.m = 0.48 + 1.574 * self.omega - 0.176 * self.omega**2 + self.a_c = 0.42748 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.08664 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + nc = self.nc + x = self.x + + # a_i(T) and da_i/dT + sqrt_Tr = np.sqrt(T / self.Tc) + alpha = (1 + self.m * (1 - sqrt_Tr))**2 + a_i = self.a_c * alpha + + # da_i/dT = a_c * 2*(1+m*(1-sqrt(T/Tc))) * m * (-1/(2*sqrt(Tc*T))) + # = -a_c * m * (1+m*(1-sqrt_Tr)) / sqrt(Tc * T) + da_i_dT = -self.a_c * self.m * (1 + self.m * (1 - sqrt_Tr)) / np.sqrt(self.Tc * T) + + # mixing rules + if nc == 1: + a_m = a_i[0] + b_m = self.b_i[0] + da_m_dT = da_i_dT[0] + else: + a_m = 0.0 + da_m_dT = 0.0 + for i in range(nc): + for j in range(nc): + a_ij = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + a_m += x[i] * x[j] * a_ij + # d(a_ij)/dT using product rule on sqrt(a_i * a_j) + if a_i[i] * a_i[j] > 0: + da_ij_dT = (1 - self.k[i, j]) / (2 * np.sqrt(a_i[i] * a_i[j])) * ( + a_i[i] * da_i_dT[j] + a_i[j] * da_i_dT[i]) + else: + da_ij_dT = 0.0 + da_m_dT += x[i] * x[j] * da_ij_dT + b_m = np.dot(x, self.b_i) + + # solve SRK cubic + A_m = a_m * P / (R**2 * T**2) + B_m = b_m * P / (R * T) + coeffs = [1, -1, A_m - B_m - B_m**2, -A_m * B_m] + Z = _solve_cubic_eos(coeffs, self.phase) + + # enthalpy departure + dh = R * T * (Z - 1) + (T * da_m_dT - a_m) / b_m * np.log(1 + b_m * P / (Z * R * T)) + return dh + + +class EnthalpyDeparturePR(Function): + r"""Isothermal enthalpy departure from the Peng-Robinson EoS (IK-CAPE Chapter 9.7.2). + + Computes the departure of the real-gas enthalpy from the ideal-gas + enthalpy using the Peng-Robinson equation of state: + + .. math:: + + \Delta h_i = RT(Z_i - 1) - \frac{1}{2\sqrt{2}\,b_i} + \left(a_i(T) - T\frac{\partial a_i}{\partial T}\right) + \ln\frac{v_i + (1+\sqrt{2})b_i}{v_i + (1-\sqrt{2})b_i} + + For a mixture, subscript *i* is replaced by mixture quantities + :math:`a_m, b_m, Z_m, v_m`. + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output port:** ``dh`` -- enthalpy departure [J/mol]. + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. + x : array_like, optional + Mole fractions [N]. Required for mixtures. + k : array_like, optional + Binary interaction parameters [N x N]. Default: zeros. + phase : str, optional + ``"vapor"`` (default) or ``"liquid"``. + """ + + input_port_labels = {"T": 0, "P": 1} + output_port_labels = {"dh": 0} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.nc = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.nc == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.nc, self.nc)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.nc, self.nc) + + # PR constants + self.m = 0.37464 + 1.54226 * self.omega - 0.26992 * self.omega**2 + self.a_c = 0.45724 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.07780 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + nc = self.nc + x = self.x + sqrt2 = np.sqrt(2.0) + + # a_i(T) and da_i/dT + sqrt_Tr = np.sqrt(T / self.Tc) + alpha = (1 + self.m * (1 - sqrt_Tr))**2 + a_i = self.a_c * alpha + da_i_dT = -self.a_c * self.m * (1 + self.m * (1 - sqrt_Tr)) / np.sqrt(self.Tc * T) + + # mixing rules + if nc == 1: + a_m = a_i[0] + b_m = self.b_i[0] + da_m_dT = da_i_dT[0] + else: + a_m = 0.0 + da_m_dT = 0.0 + for i in range(nc): + for j in range(nc): + a_ij = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + a_m += x[i] * x[j] * a_ij + if a_i[i] * a_i[j] > 0: + da_ij_dT = (1 - self.k[i, j]) / (2 * np.sqrt(a_i[i] * a_i[j])) * ( + a_i[i] * da_i_dT[j] + a_i[j] * da_i_dT[i]) + else: + da_ij_dT = 0.0 + da_m_dT += x[i] * x[j] * da_ij_dT + b_m = np.dot(x, self.b_i) + + # solve PR cubic + A_m = a_m * P / (R**2 * T**2) + B_m = b_m * P / (R * T) + coeffs = [1, + -(1 - B_m), + A_m - 3 * B_m**2 - 2 * B_m, + -(A_m * B_m - B_m**2 - B_m**3)] + Z = _solve_cubic_eos(coeffs, self.phase) + v = Z * R * T / P + + # enthalpy departure (PR) + log_term = np.log((v + (1 + sqrt2) * b_m) / (v + (1 - sqrt2) * b_m)) + dh = R * T * (Z - 1) - 1.0 / (2 * sqrt2 * b_m) * (a_m - T * da_m_dT) * log_term + return dh diff --git a/src/pathsim_chem/thermodynamics/fugacity_coefficients.py b/src/pathsim_chem/thermodynamics/fugacity_coefficients.py new file mode 100644 index 0000000..54e912c --- /dev/null +++ b/src/pathsim_chem/thermodynamics/fugacity_coefficients.py @@ -0,0 +1,353 @@ +######################################################################################### +## +## IK-CAPE Fugacity Coefficients (Chapter 8) +## +## Fugacity coefficients from cubic equations of state and the +## virial equation for vapor-liquid equilibrium calculations. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + +from .equations_of_state import _solve_cubic_eos + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# BLOCKS ================================================================================ + +class FugacityRKS(Function): + r"""Fugacity coefficients from the Soave-Redlich-Kwong EoS (IK-CAPE Chapter 8.1). + + Computes the fugacity coefficient of each component in a vapor or liquid + phase using the Soave-Redlich-Kwong cubic equation of state. The fugacity + coefficient :math:`\phi_i` relates the fugacity of a real gas to its + partial pressure and is essential for rigorous VLE calculations: + + .. math:: + + f_i = \phi_i \, y_i \, P + + For a pure component the fugacity coefficient is: + + .. math:: + + \ln \phi_i = Z_i - 1 - \ln(Z_i - B_i) + - \frac{A_i}{B_i} \ln\!\left(1 + \frac{B_i}{Z_i}\right) + + For component *i* in a mixture: + + .. math:: + + \ln \phi_i = \frac{b_i}{b_m}(Z_m - 1) - \ln(Z_m - B_m) + + \frac{A_m}{B_m}\left(\frac{b_i}{b_m} + - \frac{2\sum_j x_j a_{ij}}{a_m}\right) + \ln\!\left(1 + \frac{B_m}{Z_m}\right) + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output ports:** ``phi_0``, ``phi_1``, ... -- fugacity coefficient [-] + for each component (one output per component). + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. + x : array_like, optional + Mole fractions [N]. Required for mixtures, omit for pure components. + k : array_like, optional + Binary interaction parameters [N x N]. Default: zero for all pairs. + phase : str, optional + ``"vapor"`` (default) or ``"liquid"`` -- selects the EoS root. + """ + + input_port_labels = {"T": 0, "P": 1} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.nc = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.nc == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.nc, self.nc)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.nc, self.nc) + + # constant pure component parameters (SRK) + self.m = 0.48 + 1.574 * self.omega - 0.176 * self.omega**2 + self.a_c = 0.42748 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.08664 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + nc = self.nc + x = self.x + + # temperature-dependent a_i + alpha = (1 + self.m * (1 - np.sqrt(T / self.Tc)))**2 + a_i = self.a_c * alpha + + # cross parameters a_ij + a_ij = np.zeros((nc, nc)) + for i in range(nc): + for j in range(nc): + a_ij[i, j] = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + + # mixing rules + if nc == 1: + a_m = a_i[0] + b_m = self.b_i[0] + else: + a_m = 0.0 + for i in range(nc): + for j in range(nc): + a_m += x[i] * x[j] * a_ij[i, j] + b_m = np.dot(x, self.b_i) + + # dimensionless parameters + A_m = a_m * P / (R**2 * T**2) + B_m = b_m * P / (R * T) + + # solve cubic: Z^3 - Z^2 + (A-B-B^2)*Z - AB = 0 + coeffs = [1, -1, A_m - B_m - B_m**2, -A_m * B_m] + Z = _solve_cubic_eos(coeffs, self.phase) + + # fugacity coefficients for each component + ln_phi = np.zeros(nc) + for i in range(nc): + # sum_j x_j * a_ij + sum_xa = np.dot(x, a_ij[i, :]) + ln_phi[i] = (self.b_i[i] / b_m * (Z - 1) + - np.log(Z - B_m) + + A_m / B_m * (self.b_i[i] / b_m + - 2.0 * sum_xa / a_m) + * np.log(1 + B_m / Z)) + + return tuple(np.exp(ln_phi)) + + +class FugacityPR(Function): + r"""Fugacity coefficients from the Peng-Robinson EoS (IK-CAPE Chapter 8.2). + + Computes the fugacity coefficient of each component in a vapor or liquid + phase using the Peng-Robinson cubic equation of state. The PR EoS + provides improved liquid density predictions compared to SRK and is + widely used for hydrocarbon and industrial fluid systems. + + For a pure component: + + .. math:: + + \ln \phi_i = Z_i - 1 - \ln(Z_i - B_i) + - \frac{A_i}{2\sqrt{2}\,B_i} + \ln\frac{Z_i + (1+\sqrt{2})B_i}{Z_i + (1-\sqrt{2})B_i} + + For component *i* in a mixture: + + .. math:: + + \ln \phi_i = \frac{b_i}{b_m}(Z_m - 1) - \ln(Z_m - B_m) + - \frac{A_m}{2\sqrt{2}\,B_m} + \left(\frac{2\sum_j x_j a_{ij}}{a_m} - \frac{b_i}{b_m}\right) + \ln\frac{Z_m + (1+\sqrt{2})B_m}{Z_m + (1-\sqrt{2})B_m} + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output ports:** ``phi_0``, ``phi_1``, ... -- fugacity coefficient [-] + for each component. + + Parameters + ---------- + Tc : float or array_like + Critical temperature(s) [K]. + Pc : float or array_like + Critical pressure(s) [Pa]. + omega : float or array_like + Acentric factor(s) [-]. + x : array_like, optional + Mole fractions [N]. Required for mixtures, omit for pure components. + k : array_like, optional + Binary interaction parameters [N x N]. Default: zero for all pairs. + phase : str, optional + ``"vapor"`` (default) or ``"liquid"`` -- selects the EoS root. + """ + + input_port_labels = {"T": 0, "P": 1} + + def __init__(self, Tc, Pc, omega, x=None, k=None, phase="vapor"): + self.Tc = np.atleast_1d(np.asarray(Tc, dtype=float)) + self.Pc = np.atleast_1d(np.asarray(Pc, dtype=float)) + self.omega = np.atleast_1d(np.asarray(omega, dtype=float)) + self.nc = len(self.Tc) + self.phase = phase + + if x is None: + self.x = np.ones(1) if self.nc == 1 else None + else: + self.x = np.asarray(x, dtype=float) + + if k is None: + self.k = np.zeros((self.nc, self.nc)) + else: + self.k = np.asarray(k, dtype=float).reshape(self.nc, self.nc) + + # constant pure component parameters (PR) + self.m = 0.37464 + 1.54226 * self.omega - 0.26992 * self.omega**2 + self.a_c = 0.45724 * R**2 * self.Tc**2 / self.Pc + self.b_i = 0.07780 * R * self.Tc / self.Pc + + super().__init__(func=self._eval) + + def _eval(self, T, P): + nc = self.nc + x = self.x + sqrt2 = np.sqrt(2.0) + + # temperature-dependent a_i + alpha = (1 + self.m * (1 - np.sqrt(T / self.Tc)))**2 + a_i = self.a_c * alpha + + # cross parameters a_ij + a_ij = np.zeros((nc, nc)) + for i in range(nc): + for j in range(nc): + a_ij[i, j] = np.sqrt(a_i[i] * a_i[j]) * (1 - self.k[i, j]) + + # mixing rules + if nc == 1: + a_m = a_i[0] + b_m = self.b_i[0] + else: + a_m = 0.0 + for i in range(nc): + for j in range(nc): + a_m += x[i] * x[j] * a_ij[i, j] + b_m = np.dot(x, self.b_i) + + # dimensionless parameters + A_m = a_m * P / (R**2 * T**2) + B_m = b_m * P / (R * T) + + # solve PR cubic: Z^3 - (1-B)Z^2 + (A-3B^2-2B)Z - (AB-B^2-B^3) = 0 + coeffs = [1, + -(1 - B_m), + A_m - 3 * B_m**2 - 2 * B_m, + -(A_m * B_m - B_m**2 - B_m**3)] + Z = _solve_cubic_eos(coeffs, self.phase) + + # fugacity coefficients for each component + ln_phi = np.zeros(nc) + log_term = np.log((Z + (1 + sqrt2) * B_m) / (Z + (1 - sqrt2) * B_m)) + + for i in range(nc): + sum_xa = np.dot(x, a_ij[i, :]) + ln_phi[i] = (self.b_i[i] / b_m * (Z - 1) + - np.log(Z - B_m) + - A_m / (2 * sqrt2 * B_m) + * (2.0 * sum_xa / a_m - self.b_i[i] / b_m) + * log_term) + + return tuple(np.exp(ln_phi)) + + +class FugacityVirial(Function): + r"""Fugacity coefficients from the truncated virial equation (IK-CAPE Chapter 8.3). + + Computes fugacity coefficients using the second virial equation of state, + which is valid at low to moderate pressures (typically below 10-15 bar). + The virial EoS is theoretically exact in the limit of low density and is + particularly useful for gases at moderate conditions where cubic EoS may + be unnecessarily complex. + + For a pure component: + + .. math:: + + \phi_i^0 = \frac{\exp(2 B_{ii} / v)}{1 + B_{ii}/v} + + where :math:`v` is the molar volume from: + + .. math:: + + v = \frac{RT}{2P}\left(1 + \sqrt{1 + \frac{4 P B_{ii}}{RT}}\right) + + For component *i* in a mixture: + + .. math:: + + \phi_i^v = \frac{\exp\!\left(\frac{2}{v}\sum_j y_j B_{ij}\right)} + {1 + B_m/v} + + where :math:`B_m = \sum_i \sum_j y_i y_j B_{ij}` and + :math:`Z_m = 1 + B_m/v`. + + **Input ports:** ``T`` -- temperature [K], ``P`` -- pressure [Pa]. + + **Output ports:** ``phi_0``, ``phi_1``, ... -- fugacity coefficient [-] + for each component. + + Parameters + ---------- + B : array_like + Second virial coefficients [N x N] in [m^3/mol]. For a pure component, + pass a scalar or 1x1 array. For a mixture, pass the full symmetric + matrix including cross-coefficients :math:`B_{ij}`. + These are typically functions of temperature, but here assumed constant + at the evaluation temperature. + y : array_like, optional + Vapor-phase mole fractions [N]. Required for mixtures. + """ + + input_port_labels = {"T": 0, "P": 1} + + def __init__(self, B, y=None): + self.B = np.atleast_2d(np.asarray(B, dtype=float)) + self.nc = self.B.shape[0] + + if y is None: + self.y = np.ones(1) if self.nc == 1 else None + else: + self.y = np.asarray(y, dtype=float) + + super().__init__(func=self._eval) + + def _eval(self, T, P): + nc = self.nc + y = self.y + B_mat = self.B + + # mixture second virial coefficient + B_m = 0.0 + for i in range(nc): + for j in range(nc): + B_m += y[i] * y[j] * B_mat[i, j] + + # molar volume from truncated virial + v = R * T / (2 * P) * (1 + np.sqrt(1 + 4 * P * B_m / (R * T))) + + # fugacity coefficients + ln_phi = np.zeros(nc) + for i in range(nc): + sum_yB = np.dot(y, B_mat[i, :]) + ln_phi[i] = 2.0 * sum_yB / v - np.log(1 + B_m / v) + + return tuple(np.exp(ln_phi)) diff --git a/src/pathsim_chem/thermodynamics/reactions.py b/src/pathsim_chem/thermodynamics/reactions.py new file mode 100644 index 0000000..480a892 --- /dev/null +++ b/src/pathsim_chem/thermodynamics/reactions.py @@ -0,0 +1,228 @@ +######################################################################################### +## +## IK-CAPE Chemical Reactions (Chapter 10) +## +## Temperature-dependent equilibrium constants, kinetic rate constants, +## and reaction rate expressions for chemical process simulation. +## +######################################################################################### + +# IMPORTS =============================================================================== + +import numpy as np + +from pathsim.blocks.function import Function + + +# CONSTANTS ============================================================================= + +R = 8.314462 # gas constant [J/(mol*K)] + + +# BLOCKS ================================================================================ + +class EquilibriumConstant(Function): + r"""Temperature-dependent chemical equilibrium constant (IK-CAPE Chapter 10.1). + + Computes the equilibrium constant :math:`K_{eq}(T)` for a chemical + reaction as a function of temperature using a six-parameter correlation. + The equilibrium constant is used in equilibrium reactor models to + determine the extent of reaction at thermodynamic equilibrium. + + It is related to the standard Gibbs energy of reaction: + + .. math:: + + \Delta G^\circ = -RT \ln K_{eq} + + Temperature dependence: + + .. math:: + + \ln K_{eq}(T) = a_0 + \frac{a_1}{T} + a_2 \ln T + + a_3 T + a_4 T^2 + a_5 T^3 + + The simplest form (van't Hoff) uses only :math:`a_0` and :math:`a_1`, + where :math:`a_1 = -\Delta H^\circ / R`. Additional terms capture + the temperature dependence of the heat of reaction via heat capacity. + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``Keq`` -- equilibrium constant [-] (dimensionless + or in appropriate pressure/concentration units depending on convention). + + Parameters + ---------- + a0 : float + Constant term. + a1 : float, optional + Coefficient of :math:`1/T`. Related to the standard enthalpy of + reaction. Default: 0. + a2 : float, optional + Coefficient of :math:`\ln T`. Related to heat capacity change. + Default: 0. + a3 : float, optional + Coefficient of :math:`T`. Default: 0. + a4 : float, optional + Coefficient of :math:`T^2`. Default: 0. + a5 : float, optional + Coefficient of :math:`T^3`. Default: 0. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"Keq": 0} + + def __init__(self, a0, a1=0.0, a2=0.0, a3=0.0, a4=0.0, a5=0.0): + self.coeffs = (a0, a1, a2, a3, a4, a5) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3, a4, a5 = self.coeffs + ln_K = a0 + a1 / T + a2 * np.log(T) + a3 * T + a4 * T**2 + a5 * T**3 + return np.exp(ln_K) + + +class KineticRateConstant(Function): + r"""Temperature-dependent kinetic rate constant (IK-CAPE Chapter 10.2). + + Computes the reaction rate constant :math:`k(T)` for a kinetically + controlled reaction using an extended Arrhenius-type correlation. + This is the fundamental expression for reaction kinetics, describing + how fast a chemical reaction proceeds at a given temperature. + + Two forms are supported: + + **Standard form** (4 parameters): + + .. math:: + + \ln k(T) = a_0 + \frac{a_1}{T} + a_2 \ln T + a_3 T + + where :math:`a_0 = \ln A` (pre-exponential factor) and + :math:`a_1 = -E_a/R` (activation energy). + + **Extended form** with separate forward and reverse constants, + where the reverse constant can be derived from: + + .. math:: + + k_{reverse} = k_{forward} / K_{eq} + + **Input port:** ``T`` -- temperature [K]. + + **Output port:** ``k`` -- rate constant (units depend on reaction order). + + Parameters + ---------- + a0 : float + Constant term (:math:`\ln A` for Arrhenius). + a1 : float, optional + Coefficient of :math:`1/T` (:math:`-E_a/R`). Default: 0. + a2 : float, optional + Coefficient of :math:`\ln T` (power-law temperature exponent). + Default: 0. + a3 : float, optional + Coefficient of :math:`T`. Default: 0. + """ + + input_port_labels = {"T": 0} + output_port_labels = {"k": 0} + + def __init__(self, a0, a1=0.0, a2=0.0, a3=0.0): + self.coeffs = (a0, a1, a2, a3) + super().__init__(func=self._eval) + + def _eval(self, T): + a0, a1, a2, a3 = self.coeffs + ln_k = a0 + a1 / T + a2 * np.log(T) + a3 * T + return np.exp(ln_k) + + +class PowerLawRate(Function): + r"""Power-law reaction rate expression (IK-CAPE Chapter 10.2 KILM/KIVM). + + Computes the volumetric rate of reaction for a system of :math:`N_r` + reactions involving :math:`N_c` species using power-law kinetics. + Each reaction rate is the product of the rate constant and the + concentrations (or mole fractions) raised to their stoichiometric powers: + + .. math:: + + r_k = f_k \prod_i c_i^{\nu_{i,k}} + + For equilibrium-limited reactions, the rate is modified by an + equilibrium driving-force term: + + .. math:: + + r_k = f_k \prod_i c_i^{\nu_{i,k}^{fwd}} + \left(1 - \frac{\prod_i c_i^{\nu_{i,k}}}{K_{eq,k}}\right) + + This block computes the **net production rate** of each species as a + linear combination of individual reaction rates: + + .. math:: + + R_i = \sum_k \nu_{i,k} \, r_k + + **Input ports:** + + - ``f_0``, ``f_1``, ... -- rate constants for each reaction (from + :class:`KineticRateConstant` blocks upstream) + - ``c_0``, ``c_1``, ... -- concentrations or mole fractions of species + + **Output ports:** ``R_0``, ``R_1``, ... -- net production rate for each + species [mol/(m^3 s)] or [1/s] depending on the units of *f* and *c*. + + Parameters + ---------- + nu : array_like + Stoichiometric coefficient matrix [N_c x N_r]. Negative for + reactants, positive for products. + nu_fwd : array_like, optional + Forward reaction orders [N_c x N_r]. If not given, the absolute + values of the negative stoichiometric coefficients are used (i.e., + elementary reaction kinetics). + Keq : array_like, optional + Equilibrium constants [N_r]. If provided, an equilibrium driving + force term is included. Default: None (irreversible reactions). + """ + + def __init__(self, nu, nu_fwd=None, Keq=None): + self.nu = np.asarray(nu, dtype=float) + self.n_species, self.n_rxn = self.nu.shape + + if nu_fwd is None: + # elementary kinetics: forward orders = |negative stoich coeffs| + self.nu_fwd = np.where(self.nu < 0, -self.nu, 0.0) + else: + self.nu_fwd = np.asarray(nu_fwd, dtype=float) + + self.Keq = None if Keq is None else np.asarray(Keq, dtype=float) + + super().__init__(func=self._eval) + + def _eval(self, *inputs): + nr = self.n_rxn + nc = self.n_species + + # first nr inputs are rate constants, next nc are concentrations + f = np.array(inputs[:nr]) + c = np.array(inputs[nr:nr + nc]) + + # compute reaction rates + rates = np.zeros(nr) + for k in range(nr): + # forward rate + rate_fwd = f[k] * np.prod(c ** self.nu_fwd[:, k]) + + if self.Keq is not None: + # equilibrium driving force: (1 - prod(c^nu) / Keq) + Q = np.prod(c ** self.nu[:, k]) + rate_fwd *= (1 - Q / self.Keq[k]) + + rates[k] = rate_fwd + + # net production rates for each species + R = np.dot(self.nu, rates) + return tuple(R) diff --git a/tests/thermodynamics/test_activity_coefficients.py b/tests/thermodynamics/test_activity_coefficients.py index 43c687c..f6ad11b 100644 --- a/tests/thermodynamics/test_activity_coefficients.py +++ b/tests/thermodynamics/test_activity_coefficients.py @@ -16,6 +16,7 @@ NRTL, Wilson, UNIQUAC, + FloryHuggins, ) @@ -217,6 +218,76 @@ def test_three_components(self): self.assertTrue(np.isfinite(g)) +class TestFloryHuggins(unittest.TestCase): + + def test_init(self): + FH = FloryHuggins( + x=[0.8, 0.2], + r=[100.0], + chi0=[0.45], + ) + self.assertEqual(FH.n, 2) + + def test_binary_pure_solvent(self): + # Pure solvent: gamma_1 should be 1.0 + FH = FloryHuggins( + x=[1.0, 0.0], + r=[100.0], + chi0=[0.45], + ) + gamma1, gamma2 = eval_block(FH, 300) + self.assertAlmostEqual(gamma1, 1.0, places=10) + + def test_binary_solvent_polymer(self): + # Solvent-polymer system with typical polystyrene-toluene values + FH = FloryHuggins( + x=[0.9, 0.1], + r=[100.0], + chi0=[0.45], + ) + gamma1, gamma2 = eval_block(FH, 300) + self.assertGreater(gamma1, 0) + self.assertGreater(gamma2, 0) + self.assertTrue(np.isfinite(gamma1)) + self.assertTrue(np.isfinite(gamma2)) + + def test_temperature_dependence(self): + FH = FloryHuggins( + x=[0.8, 0.2], + r=[50.0], + chi0=[0.3], + chi1=[100.0], + ) + g300 = eval_block(FH, 300) + g400 = eval_block(FH, 400) + self.assertNotAlmostEqual(g300[0], g400[0], places=3) + + def test_multiple_polymers(self): + # Solvent + two polymer species + FH = FloryHuggins( + x=[0.7, 0.2, 0.1], + r=[50.0, 200.0], + chi0=[0.4, 0.6], + ) + gammas = eval_block(FH, 300) + self.assertEqual(len(gammas), 3) + for g in gammas: + self.assertGreater(g, 0) + self.assertTrue(np.isfinite(g)) + + def test_zero_chi_ideal_entropy(self): + # With chi=0, only combinatorial (entropic) contributions + FH = FloryHuggins( + x=[0.5, 0.5], + r=[10.0], + chi0=[0.0], + ) + gamma1, gamma2 = eval_block(FH, 300) + # gamma_1 should still deviate from 1 due to entropy of mixing + self.assertGreater(gamma1, 0) + self.assertTrue(np.isfinite(gamma1)) + + # SMOKE TEST =========================================================================== class TestSmokeAllActivityCoefficients(unittest.TestCase): @@ -226,6 +297,7 @@ def test_all_at_350K(self): NRTL(x=[0.5, 0.5], a=[[0, 1.0], [1.5, 0]], c=[[0, 0.3], [0.3, 0]]), Wilson(x=[0.5, 0.5], a=[[0, 0.5], [-0.3, 0]]), UNIQUAC(x=[0.5, 0.5], r=[2.1, 0.92], q=[1.97, 1.4], a=[[0, -1.0], [2.0, 0]]), + FloryHuggins(x=[0.8, 0.2], r=[100.0], chi0=[0.45]), ] for block in blocks: diff --git a/tests/thermodynamics/test_enthalpy.py b/tests/thermodynamics/test_enthalpy.py new file mode 100644 index 0000000..09c8740 --- /dev/null +++ b/tests/thermodynamics/test_enthalpy.py @@ -0,0 +1,362 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.enthalpy.py' +## +## IK-CAPE Enthalpy (Chapter 9) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + ExcessEnthalpyNRTL, + ExcessEnthalpyUNIQUAC, + ExcessEnthalpyWilson, + ExcessEnthalpyRedlichKister, + EnthalpyDepartureRKS, + EnthalpyDeparturePR, +) + + +# CONSTANTS ============================================================================ + +R = 8.314462 + + +# HELPERS ============================================================================== + +def eval_block_T(block, T): + """Set input T, update, return output.""" + block.inputs[0] = T + block.update(None) + return block.outputs[0] + + +def eval_block_TP(block, T, P): + """Set inputs T and P, update, return output.""" + block.inputs[0] = T + block.inputs[1] = P + block.update(None) + return block.outputs[0] + + +# TESTS: EXCESS ENTHALPY =============================================================== + +class TestExcessEnthalpyNRTL(unittest.TestCase): + + def test_init(self): + hE = ExcessEnthalpyNRTL( + x=[0.5, 0.5], + a=[[0, 1.5], [2.0, 0]], + ) + self.assertEqual(hE.n, 2) + + def test_zero_for_pure_component(self): + # Pure component should have hE = 0 + hE = ExcessEnthalpyNRTL( + x=[1.0, 0.0], + a=[[0, 1.5], [2.0, 0]], + b=[[0, 100], [-50, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + result = eval_block_T(hE, 350) + self.assertAlmostEqual(result, 0.0, places=5) + + def test_finite_for_mixture(self): + hE = ExcessEnthalpyNRTL( + x=[0.4, 0.6], + a=[[0, -0.801], [3.458, 0]], + b=[[0, 300], [-200, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + result = eval_block_T(hE, 350) + self.assertTrue(np.isfinite(result)) + + def test_temperature_dependence(self): + hE = ExcessEnthalpyNRTL( + x=[0.5, 0.5], + a=[[0, 1.0], [1.0, 0]], + b=[[0, 500], [500, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + h300 = eval_block_T(hE, 300) + h400 = eval_block_T(hE, 400) + self.assertNotAlmostEqual(h300, h400, places=0) + + def test_zero_when_no_temperature_terms(self): + # With only constant tau (b=0, e=0, f=0), tau'=0 + # and if d=0 (S'=0), then G'=0, so hE should be 0 + hE = ExcessEnthalpyNRTL( + x=[0.5, 0.5], + a=[[0, 1.0], [1.5, 0]], + c=[[0, 0.3], [0.3, 0]], + ) + result = eval_block_T(hE, 350) + self.assertAlmostEqual(result, 0.0, places=8) + + +class TestExcessEnthalpyUNIQUAC(unittest.TestCase): + + def test_init(self): + hE = ExcessEnthalpyUNIQUAC( + x=[0.5, 0.5], + r=[2.1055, 0.92], + q=[1.972, 1.4], + a=[[0, -1.318], [2.772, 0]], + ) + self.assertEqual(hE.n, 2) + + def test_finite_for_mixture(self): + hE = ExcessEnthalpyUNIQUAC( + x=[0.4, 0.6], + r=[2.1055, 0.92], + q=[1.972, 1.4], + a=[[0, -1.318], [2.772, 0]], + b=[[0, 200], [-100, 0]], + ) + result = eval_block_T(hE, 350) + self.assertTrue(np.isfinite(result)) + + def test_zero_when_no_temperature_terms(self): + # With only constant tau exponent (b=0, c=0, d=0), tau'=0, hE=0 + hE = ExcessEnthalpyUNIQUAC( + x=[0.5, 0.5], + r=[2.1, 0.92], + q=[1.97, 1.4], + a=[[0, -1.0], [2.0, 0]], + ) + result = eval_block_T(hE, 350) + self.assertAlmostEqual(result, 0.0, places=8) + + def test_temperature_dependence(self): + hE = ExcessEnthalpyUNIQUAC( + x=[0.5, 0.5], + r=[2.1, 0.92], + q=[1.97, 1.4], + a=[[0, -1.0], [2.0, 0]], + b=[[0, 500], [-300, 0]], + ) + h300 = eval_block_T(hE, 300) + h400 = eval_block_T(hE, 400) + self.assertNotAlmostEqual(h300, h400, places=0) + + +class TestExcessEnthalpyWilson(unittest.TestCase): + + def test_init(self): + hE = ExcessEnthalpyWilson( + x=[0.5, 0.5], + a=[[0, 0.5], [-0.3, 0]], + ) + self.assertEqual(hE.n, 2) + + def test_zero_when_no_temperature_terms(self): + # With only constant Lambda (b=0, c=0, d=0), Lambda'=0, hE=0 + hE = ExcessEnthalpyWilson( + x=[0.5, 0.5], + a=[[0, 0.5], [-0.3, 0]], + ) + result = eval_block_T(hE, 350) + self.assertAlmostEqual(result, 0.0, places=8) + + def test_finite_for_mixture(self): + hE = ExcessEnthalpyWilson( + x=[0.4, 0.6], + a=[[0, 0.5], [-0.3, 0]], + b=[[0, 200], [-150, 0]], + ) + result = eval_block_T(hE, 350) + self.assertTrue(np.isfinite(result)) + + def test_temperature_dependence(self): + hE = ExcessEnthalpyWilson( + x=[0.5, 0.5], + a=[[0, 0.5], [-0.3, 0]], + b=[[0, 500], [-300, 0]], + ) + h300 = eval_block_T(hE, 300) + h400 = eval_block_T(hE, 400) + self.assertNotAlmostEqual(h300, h400, places=0) + + def test_ideal_solution(self): + # With a=0, Lambda=1, Lambda'=0 => hE=0 + hE = ExcessEnthalpyWilson( + x=[0.5, 0.5], + a=[[0, 0], [0, 0]], + ) + result = eval_block_T(hE, 350) + self.assertAlmostEqual(result, 0.0, places=10) + + +class TestExcessEnthalpyRedlichKister(unittest.TestCase): + + def test_init(self): + hE = ExcessEnthalpyRedlichKister( + x=[0.5, 0.5], + coeffs={(0, 1): [[1000.0]]}, + ) + 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 + 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) + + 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 + hE = ExcessEnthalpyRedlichKister( + x=[0.3, 0.7], + coeffs={(0, 1): [[2000.0]]}, + ) + result = eval_block_T(hE, 300) + expected = 0.5 * 0.3 * 0.7 / 1.0 * 2000.0 * (0.3 - 0.7) + self.assertAlmostEqual(result, expected, places=5) + + def test_temperature_dependent_coeffs(self): + # A(T) = 1000 + 2*T + hE = ExcessEnthalpyRedlichKister( + x=[0.4, 0.6], + coeffs={(0, 1): [[1000.0, 2.0]]}, + ) + h300 = eval_block_T(hE, 300) + h400 = eval_block_T(hE, 400) + self.assertNotAlmostEqual(h300, h400, places=0) + + def test_finite(self): + hE = ExcessEnthalpyRedlichKister( + x=[0.4, 0.6], + coeffs={(0, 1): [[1000.0], [500.0], [200.0]]}, + ) + result = eval_block_T(hE, 350) + self.assertTrue(np.isfinite(result)) + + +# TESTS: ENTHALPY DEPARTURE ============================================================ + +class TestEnthalpyDepartureRKS(unittest.TestCase): + + def test_init(self): + dH = EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(dH.nc, 1) + + def test_ideal_gas_limit(self): + # At very low P, departure should be near 0 + dH = EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 1000, 100) + self.assertAlmostEqual(result, 0.0, delta=10) + + def test_finite(self): + dH = EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 300, 101325) + self.assertTrue(np.isfinite(result)) + + def test_negative_departure_vapor(self): + # For real gas at moderate P, departure is typically negative + dH = EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 300, 3e6) + self.assertLess(result, 0) + + def test_mixture(self): + dH = EnthalpyDepartureRKS( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + result = eval_block_TP(dH, 300, 101325) + self.assertTrue(np.isfinite(result)) + + +class TestEnthalpyDeparturePR(unittest.TestCase): + + def test_init(self): + dH = EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(dH.nc, 1) + + def test_ideal_gas_limit(self): + dH = EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 1000, 100) + self.assertAlmostEqual(result, 0.0, delta=10) + + def test_finite(self): + dH = EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 300, 101325) + self.assertTrue(np.isfinite(result)) + + def test_negative_departure_vapor(self): + dH = EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011) + result = eval_block_TP(dH, 300, 3e6) + self.assertLess(result, 0) + + def test_pr_vs_rks_similar(self): + # At moderate conditions, RKS and PR should give similar departure + dH_pr = EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011) + dH_rks = EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + h_pr = eval_block_TP(dH_pr, 300, 101325) + h_rks = eval_block_TP(dH_rks, 300, 101325) + # Within 20% of each other + self.assertAlmostEqual(h_pr, h_rks, delta=abs(h_rks) * 0.2 + 1) + + def test_mixture(self): + dH = EnthalpyDeparturePR( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + result = eval_block_TP(dH, 300, 101325) + self.assertTrue(np.isfinite(result)) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllEnthalpy(unittest.TestCase): + + def test_excess_enthalpy_blocks(self): + blocks = [ + ExcessEnthalpyNRTL( + x=[0.5, 0.5], a=[[0, 1.0], [1.5, 0]], + b=[[0, 100], [100, 0]], c=[[0, 0.3], [0.3, 0]]), + ExcessEnthalpyUNIQUAC( + x=[0.5, 0.5], r=[2.1, 0.92], q=[1.97, 1.4], + a=[[0, -1.0], [2.0, 0]], b=[[0, 200], [-100, 0]]), + ExcessEnthalpyWilson( + x=[0.5, 0.5], a=[[0, 0.5], [-0.3, 0]], + b=[[0, 200], [-150, 0]]), + ExcessEnthalpyRedlichKister( + x=[0.4, 0.6], coeffs={(0, 1): [[1000.0, 2.0]]}), + ] + + for block in blocks: + with self.subTest(block=block.__class__.__name__): + result = eval_block_T(block, 350) + self.assertTrue(np.isfinite(result), + f"{block.__class__.__name__} non-finite") + + def test_departure_blocks(self): + blocks = [ + EnthalpyDepartureRKS(Tc=190.6, Pc=4.6e6, omega=0.011), + EnthalpyDeparturePR(Tc=190.6, Pc=4.6e6, omega=0.011), + ] + + for block in blocks: + with self.subTest(block=block.__class__.__name__): + result = eval_block_TP(block, 300, 101325) + self.assertTrue(np.isfinite(result)) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) diff --git a/tests/thermodynamics/test_fugacity_coefficients.py b/tests/thermodynamics/test_fugacity_coefficients.py new file mode 100644 index 0000000..4e71bea --- /dev/null +++ b/tests/thermodynamics/test_fugacity_coefficients.py @@ -0,0 +1,195 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.fugacity_coefficients.py' +## +## IK-CAPE Fugacity Coefficients (Chapter 8) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + FugacityRKS, + FugacityPR, + FugacityVirial, +) + + +# CONSTANTS ============================================================================ + +R = 8.314462 + + +# HELPERS ============================================================================== + +def eval_block(block, T, P): + """Set T and P inputs, update, return outputs.""" + block.inputs[0] = T + block.inputs[1] = P + block.update(None) + n_out = len([k for k in block.outputs]) + if n_out == 1: + return (block.outputs[0],) + return tuple(block.outputs[i] for i in range(n_out)) + + +# TESTS ================================================================================ + +class TestFugacityRKS(unittest.TestCase): + + def test_init_pure(self): + F = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(F.nc, 1) + + def test_ideal_gas_limit(self): + # At very low P and high T, phi -> 1 + F = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + phi = eval_block(F, 1000, 100) + self.assertAlmostEqual(phi[0], 1.0, delta=0.01) + + def test_methane_moderate(self): + # Methane at 300K, 1 atm - phi close to 1 + F = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + phi = eval_block(F, 300, 101325) + self.assertAlmostEqual(phi[0], 1.0, delta=0.05) + self.assertGreater(phi[0], 0) + + def test_high_pressure(self): + # At higher pressure, phi deviates from 1 + F = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + phi_low = eval_block(F, 300, 101325) + phi_high = eval_block(F, 300, 3e6) + # phi should deviate more at high pressure + self.assertGreater(abs(phi_high[0] - 1.0), abs(phi_low[0] - 1.0)) + + def test_mixture(self): + # Methane + ethane mixture + F = FugacityRKS( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + phis = eval_block(F, 300, 101325) + self.assertEqual(len(phis), 2) + for phi in phis: + self.assertGreater(phi, 0) + self.assertTrue(np.isfinite(phi)) + + def test_liquid_phase(self): + F = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid") + phi = eval_block(F, 120, 5e6) + self.assertGreater(phi[0], 0) + self.assertTrue(np.isfinite(phi[0])) + + +class TestFugacityPR(unittest.TestCase): + + def test_init_pure(self): + F = FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011) + self.assertEqual(F.nc, 1) + + def test_ideal_gas_limit(self): + F = FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011) + phi = eval_block(F, 1000, 100) + self.assertAlmostEqual(phi[0], 1.0, delta=0.01) + + def test_methane_moderate(self): + F = FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011) + phi = eval_block(F, 300, 101325) + self.assertAlmostEqual(phi[0], 1.0, delta=0.05) + + def test_mixture(self): + F = FugacityPR( + Tc=[190.6, 305.3], + Pc=[4.6e6, 4.872e6], + omega=[0.011, 0.099], + x=[0.7, 0.3], + ) + phis = eval_block(F, 300, 101325) + self.assertEqual(len(phis), 2) + for phi in phis: + self.assertGreater(phi, 0) + self.assertTrue(np.isfinite(phi)) + + def test_pr_vs_rks_similar(self): + # At ideal gas conditions, both should give similar phi + F_pr = FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011) + F_rks = FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011) + phi_pr = eval_block(F_pr, 300, 101325) + phi_rks = eval_block(F_rks, 300, 101325) + self.assertAlmostEqual(phi_pr[0], phi_rks[0], delta=0.02) + + def test_liquid_phase(self): + F = FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011, phase="liquid") + phi = eval_block(F, 120, 5e6) + self.assertGreater(phi[0], 0) + self.assertTrue(np.isfinite(phi[0])) + + +class TestFugacityVirial(unittest.TestCase): + + def test_init_pure(self): + F = FugacityVirial(B=[[-4.3e-4]]) + self.assertEqual(F.nc, 1) + + def test_ideal_gas_limit(self): + # With B=0, phi should be 1 + F = FugacityVirial(B=[[0.0]], y=[1.0]) + phi = eval_block(F, 300, 101325) + self.assertAlmostEqual(phi[0], 1.0, places=5) + + def test_negative_B(self): + # With negative B (attractive interactions), phi < 1 + F = FugacityVirial(B=[[-5e-4]], y=[1.0]) + phi = eval_block(F, 300, 101325) + self.assertLess(phi[0], 1.0) + self.assertGreater(phi[0], 0) + + def test_mixture(self): + # Binary mixture with cross virial coefficients + F = FugacityVirial( + B=[[-4.3e-4, -3.5e-4], + [-3.5e-4, -6.0e-4]], + y=[0.7, 0.3], + ) + phis = eval_block(F, 300, 101325) + self.assertEqual(len(phis), 2) + for phi in phis: + self.assertGreater(phi, 0) + self.assertTrue(np.isfinite(phi)) + + def test_low_pressure(self): + # At very low pressure, phi -> 1 even with non-zero B + F = FugacityVirial(B=[[-5e-4]], y=[1.0]) + phi = eval_block(F, 500, 10) # 10 Pa + self.assertAlmostEqual(phi[0], 1.0, delta=0.001) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllFugacity(unittest.TestCase): + + def test_all_blocks(self): + blocks = [ + FugacityRKS(Tc=190.6, Pc=4.6e6, omega=0.011), + FugacityPR(Tc=190.6, Pc=4.6e6, omega=0.011), + FugacityVirial(B=[[-4.3e-4]], y=[1.0]), + ] + + for block in blocks: + with self.subTest(block=block.__class__.__name__): + phis = eval_block(block, 300, 101325) + for phi in phis: + self.assertTrue(np.isfinite(phi)) + self.assertGreater(phi, 0) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) diff --git a/tests/thermodynamics/test_reactions.py b/tests/thermodynamics/test_reactions.py new file mode 100644 index 0000000..5d7624e --- /dev/null +++ b/tests/thermodynamics/test_reactions.py @@ -0,0 +1,224 @@ +######################################################################################## +## +## TESTS FOR +## 'thermodynamics.reactions.py' +## +## IK-CAPE Chemical Reactions (Chapter 10) +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.thermodynamics import ( + EquilibriumConstant, + KineticRateConstant, + PowerLawRate, +) + + +# CONSTANTS ============================================================================ + +R = 8.314462 + + +# HELPERS ============================================================================== + +def eval_block_T(block, T): + """Set input T, update, return output.""" + block.inputs[0] = T + block.update(None) + return block.outputs[0] + + +def eval_block_multi(block, *values): + """Set multiple inputs, update, return all outputs.""" + for i, v in enumerate(values): + block.inputs[i] = v + block.update(None) + n_out = len([k for k in block.outputs]) + if n_out == 1: + return (block.outputs[0],) + return tuple(block.outputs[i] for i in range(n_out)) + + +# TESTS ================================================================================ + +class TestEquilibriumConstant(unittest.TestCase): + + def test_init(self): + K = EquilibriumConstant(a0=10.0, a1=-5000.0) + self.assertEqual(K.coeffs, (10.0, -5000.0, 0.0, 0.0, 0.0, 0.0)) + + def test_constant_only(self): + # K = exp(a0) + K = EquilibriumConstant(a0=5.0) + result = eval_block_T(K, 300) + self.assertAlmostEqual(result, np.exp(5.0), places=5) + + def test_vant_hoff(self): + # ln K = a0 + a1/T => K = exp(a0 + a1/T) + a0, a1 = 20.0, -5000.0 + K = EquilibriumConstant(a0=a0, a1=a1) + T = 400 + expected = np.exp(a0 + a1 / T) + result = eval_block_T(K, T) + self.assertAlmostEqual(result, expected, places=5) + + def test_temperature_dependence(self): + K = EquilibriumConstant(a0=10.0, a1=-5000.0) + K_300 = eval_block_T(K, 300) + K_500 = eval_block_T(K, 500) + # With negative a1, ln K increases with T (less negative 1/T term) + self.assertGreater(K_500, K_300) + + def test_all_parameters(self): + a0, a1, a2, a3, a4, a5 = 10.0, -2000.0, 1.5, -0.001, 1e-6, -1e-9 + T = 350 + K = EquilibriumConstant(a0=a0, a1=a1, a2=a2, a3=a3, a4=a4, a5=a5) + expected = np.exp(a0 + a1/T + a2*np.log(T) + a3*T + a4*T**2 + a5*T**3) + result = eval_block_T(K, T) + self.assertAlmostEqual(result, expected, places=5) + + def test_positive(self): + K = EquilibriumConstant(a0=10.0, a1=-5000.0) + result = eval_block_T(K, 350) + self.assertGreater(result, 0) + + +class TestKineticRateConstant(unittest.TestCase): + + def test_init(self): + k = KineticRateConstant(a0=30.0, a1=-10000.0) + self.assertEqual(k.coeffs, (30.0, -10000.0, 0.0, 0.0)) + + def test_arrhenius(self): + # k = A * exp(-Ea/RT) => ln k = ln(A) - Ea/(RT) = a0 + a1/T + # where a0 = ln(A), a1 = -Ea/R + a0 = np.log(1e13) # pre-exponential factor A = 1e13 + a1 = -15000.0 # Ea/R = 15000 K + k = KineticRateConstant(a0=a0, a1=a1) + T = 500 + expected = np.exp(a0 + a1 / T) + result = eval_block_T(k, T) + self.assertAlmostEqual(result, expected, places=5) + + def test_temperature_dependence(self): + k = KineticRateConstant(a0=30.0, a1=-10000.0) + k_300 = eval_block_T(k, 300) + k_500 = eval_block_T(k, 500) + # Rate constant increases with temperature + self.assertGreater(k_500, k_300) + + def test_positive(self): + k = KineticRateConstant(a0=30.0, a1=-10000.0) + result = eval_block_T(k, 350) + self.assertGreater(result, 0) + + def test_all_parameters(self): + a0, a1, a2, a3 = 25.0, -8000.0, 0.5, -0.001 + T = 400 + k = KineticRateConstant(a0=a0, a1=a1, a2=a2, a3=a3) + expected = np.exp(a0 + a1/T + a2*np.log(T) + a3*T) + result = eval_block_T(k, T) + self.assertAlmostEqual(result, expected, places=5) + + +class TestPowerLawRate(unittest.TestCase): + + def test_simple_irreversible(self): + # A -> B, rate = k * c_A + # nu = [[-1], [1]], forward order: |nu_A| = 1 + rate = PowerLawRate(nu=[[-1], [1]]) + + # Inputs: f_0=k=2.0, c_0=c_A=0.5, c_1=c_B=0.3 + R_vals = eval_block_multi(rate, 2.0, 0.5, 0.3) + # r = k * c_A = 2.0 * 0.5 = 1.0 + # R_A = -1 * 1.0 = -1.0, R_B = 1 * 1.0 = 1.0 + self.assertAlmostEqual(R_vals[0], -1.0, places=10) + self.assertAlmostEqual(R_vals[1], 1.0, places=10) + + def test_bimolecular(self): + # A + B -> C, rate = k * c_A * c_B + rate = PowerLawRate(nu=[[-1], [-1], [1]]) + + # k=3.0, c_A=0.4, c_B=0.5, c_C=0.1 + R_vals = eval_block_multi(rate, 3.0, 0.4, 0.5, 0.1) + # r = 3.0 * 0.4 * 0.5 = 0.6 + self.assertAlmostEqual(R_vals[0], -0.6, places=10) + self.assertAlmostEqual(R_vals[1], -0.6, places=10) + self.assertAlmostEqual(R_vals[2], 0.6, places=10) + + def test_equilibrium_limited(self): + # A <-> B, Keq = 10.0 + rate = PowerLawRate( + nu=[[-1], [1]], + Keq=[10.0], + ) + + # At c_A=1.0, c_B=0.0: driving force = 1 - (c_B/c_A) / Keq = 1 + R_vals = eval_block_multi(rate, 2.0, 1.0, 0.0) + # r = k * c_A * (1 - c_B/(c_A*Keq)) = 2.0 * 1.0 * (1 - 0) = 2.0 + self.assertAlmostEqual(R_vals[0], -2.0, places=10) + + def test_equilibrium_at_equilibrium(self): + # At equilibrium, net rate should be zero + # A <-> B, Keq = 4.0, so at equilibrium c_B/c_A = 4.0 + rate = PowerLawRate( + nu=[[-1], [1]], + Keq=[4.0], + ) + # c_A=0.2, c_B=0.8 => c_B/c_A = 4.0 = Keq + R_vals = eval_block_multi(rate, 1.0, 0.2, 0.8) + self.assertAlmostEqual(R_vals[0], 0.0, places=10) + self.assertAlmostEqual(R_vals[1], 0.0, places=10) + + def test_two_reactions(self): + # A -> B (rxn 1), A -> C (rxn 2) + # nu = [[-1, -1], [1, 0], [0, 1]] + rate = PowerLawRate(nu=[[-1, -1], [1, 0], [0, 1]]) + + # k1=1.0, k2=2.0, c_A=0.5, c_B=0.2, c_C=0.3 + R_vals = eval_block_multi(rate, 1.0, 2.0, 0.5, 0.2, 0.3) + # r1 = 1.0 * 0.5 = 0.5, r2 = 2.0 * 0.5 = 1.0 + # R_A = -0.5 - 1.0 = -1.5, R_B = 0.5, R_C = 1.0 + self.assertAlmostEqual(R_vals[0], -1.5, places=10) + self.assertAlmostEqual(R_vals[1], 0.5, places=10) + self.assertAlmostEqual(R_vals[2], 1.0, places=10) + + def test_custom_orders(self): + # A -> B, but with fractional order 0.5 for A + rate = PowerLawRate( + nu=[[-1], [1]], + nu_fwd=[[0.5], [0.0]], + ) + # k=4.0, c_A=0.25, c_B=0.5 + R_vals = eval_block_multi(rate, 4.0, 0.25, 0.5) + # r = 4.0 * 0.25^0.5 = 4.0 * 0.5 = 2.0 + self.assertAlmostEqual(R_vals[0], -2.0, places=10) + self.assertAlmostEqual(R_vals[1], 2.0, places=10) + + +# SMOKE TEST =========================================================================== + +class TestSmokeAllReactions(unittest.TestCase): + + def test_all_correlation_blocks(self): + blocks_T = [ + EquilibriumConstant(a0=10.0, a1=-5000.0), + KineticRateConstant(a0=30.0, a1=-10000.0), + ] + + for block in blocks_T: + with self.subTest(block=block.__class__.__name__): + result = eval_block_T(block, 350) + self.assertTrue(np.isfinite(result)) + self.assertGreater(result, 0) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From da7ec4a1c5af2bc8ea05249ffd806d663274db1e Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 15:13:11 +0100 Subject: [PATCH 6/9] Update README and add example notebooks for docs --- README.md | 96 ++++++- docs/source/examples/equation_of_state.ipynb | 226 +++++++++++++++++ .../examples/vapor_pressure_curves.ipynb | 195 ++++++++++++++ docs/source/examples/vle_calculation.ipynb | 238 ++++++++++++++++++ 4 files changed, 752 insertions(+), 3 deletions(-) create mode 100644 docs/source/examples/equation_of_state.ipynb create mode 100644 docs/source/examples/vapor_pressure_curves.ipynb create mode 100644 docs/source/examples/vle_calculation.ipynb diff --git a/README.md b/README.md index 3d7b433..1866de8 100644 --- a/README.md +++ b/README.md @@ -3,8 +3,98 @@ PathSim-Chem Logo

------------- +

+ Chemical engineering blocks for PathSim +

+ +

+ PyPI + License +

+ +

+ Documentation • + PathSim Homepage • + GitHub +

+ +--- + +PathSim-Chem extends the [PathSim](https://github.com/pathsim/pathsim) simulation framework with blocks for chemical engineering and thermodynamic property calculations. All blocks follow the standard PathSim block interface and can be connected into simulation diagrams. + +## Features + +- **IK-CAPE Thermodynamics** — 50+ blocks implementing the DECHEMA IK-CAPE standard for thermodynamic property calculations +- **Pure Component Correlations** — Antoine, Wagner, Kirchhoff, Rackett, Aly-Lee, DIPPR, and 10 more temperature-dependent property correlations +- **Mixing Rules** — Linear, quadratic, Lorentz-Berthelot, and other calculation-of-averages rules for mixture properties +- **Activity Coefficients** — NRTL, Wilson, UNIQUAC, and Flory-Huggins models for liquid-phase non-ideality +- **Equations of State** — Peng-Robinson and Soave-Redlich-Kwong cubic EoS with mixture support +- **Fugacity Coefficients** — EoS-based and virial equation fugacity calculations +- **Excess Enthalpy & Departure** — NRTL, UNIQUAC, Wilson, Redlich-Kister excess enthalpy and EoS departure functions +- **Chemical Reactions** — Equilibrium constants, kinetic rate constants, and power-law rate expressions +- **Tritium Processing** — GLC columns, TCAP cascades, bubblers, and splitters for tritium separation + +## Install + +```bash +pip install pathsim-chem +``` + +## Quick Example + +Compute the vapor pressure of water at 100 °C using the Antoine equation: + +```python +from pathsim_chem.thermodynamics import Antoine + +# Antoine coefficients for water (NIST) +antoine = Antoine(a0=23.2256, a1=3835.18, a2=-45.343) + +# Evaluate at 373.15 K +antoine.inputs[0] = 373.15 +antoine.update(None) +P_sat = antoine.outputs[0] # ≈ 101325 Pa +``` + +Use activity coefficients in a simulation: + +```python +from pathsim_chem.thermodynamics import NRTL + +# Ethanol-water NRTL model +nrtl = NRTL( + x=[0.4, 0.6], + a=[[0, -0.801], [3.458, 0]], + c=[[0, 0.3], [0.3, 0]], +) + +# Evaluate at 350 K +nrtl.inputs[0] = 350 +nrtl.update(None) +gamma_ethanol = nrtl.outputs[0] +gamma_water = nrtl.outputs[1] +``` + +## Modules + +| Module | Description | +|---|---| +| `pathsim_chem.thermodynamics.correlations` | 16 pure component property correlations (Antoine, Wagner, etc.) | +| `pathsim_chem.thermodynamics.averages` | 10 mixing rules for calculation of averages | +| `pathsim_chem.thermodynamics.activity_coefficients` | NRTL, Wilson, UNIQUAC, Flory-Huggins | +| `pathsim_chem.thermodynamics.equations_of_state` | Peng-Robinson, Soave-Redlich-Kwong | +| `pathsim_chem.thermodynamics.corrections` | Poynting correction, Henry's law | +| `pathsim_chem.thermodynamics.fugacity_coefficients` | RKS, PR, and virial fugacity coefficients | +| `pathsim_chem.thermodynamics.enthalpy` | Excess enthalpy and enthalpy departure functions | +| `pathsim_chem.thermodynamics.reactions` | Equilibrium constants, rate constants, power-law rates | +| `pathsim_chem.tritium` | GLC, TCAP, bubbler, splitter blocks for tritium processing | + +## Learn More + +- [Documentation](https://docs.pathsim.org/chem) — API reference and examples +- [PathSim](https://github.com/pathsim/pathsim) — the core simulation framework +- [Contributing](https://docs.pathsim.org/pathsim/latest/contributing) — how to contribute -# PathSim Chemical Engineering Toolbox +## License -This toolbox extends the core pathsim package with domain specific component models for chemical engineering. \ No newline at end of file +MIT diff --git a/docs/source/examples/equation_of_state.ipynb b/docs/source/examples/equation_of_state.ipynb new file mode 100644 index 0000000..b8a9fd5 --- /dev/null +++ b/docs/source/examples/equation_of_state.ipynb @@ -0,0 +1,226 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Cubic Equations of State\n", + "\n", + "Computing compressibility factors and molar volumes using the Peng-Robinson and Soave-Redlich-Kwong cubic equations of state for pure methane and a methane-ethane mixture." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Import the equation of state blocks." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim_chem.thermodynamics import PengRobinson, RedlichKwongSoave" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pure Methane\n", + "\n", + "The Peng-Robinson EoS solves a cubic equation in the compressibility factor $Z = Pv/(RT)$:\n", + "\n", + "$$Z^3 - (1-B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0$$\n", + "\n", + "where $A$ and $B$ are dimensionless parameters computed from the critical properties. Each block takes $(T, P)$ as inputs and returns $(v, Z)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Critical properties of methane\n", + "Tc_CH4 = 190.6 # K\n", + "Pc_CH4 = 4.6e6 # Pa\n", + "omega_CH4 = 0.011 # acentric factor\n", + "\n", + "# Create EoS blocks\n", + "pr = PengRobinson(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", + "rks = RedlichKwongSoave(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", + "\n", + "def eval_eos(block, T, P):\n", + " \"\"\"Evaluate an EoS block and return (v, Z).\"\"\"\n", + " block.inputs[0] = T\n", + " block.inputs[1] = P\n", + " block.update(None)\n", + " return block.outputs[0], block.outputs[1]\n", + "\n", + "# Evaluate at 300 K, 1 atm\n", + "v_pr, Z_pr = eval_eos(pr, 300, 101325)\n", + "v_rks, Z_rks = eval_eos(rks, 300, 101325)\n", + "\n", + "print(f\"Methane at 300 K, 1 atm:\")\n", + "print(f\" PR: Z = {Z_pr:.6f}, v = {v_pr:.6e} m³/mol\")\n", + "print(f\" RKS: Z = {Z_rks:.6f}, v = {v_rks:.6e} m³/mol\")\n", + "print(f\" Ideal gas: Z = 1.0, v = {8.314462 * 300 / 101325:.6e} m³/mol\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compressibility Factor vs Pressure\n", + "\n", + "At low pressure, $Z \\to 1$ (ideal gas). As pressure increases, intermolecular forces cause $Z$ to deviate. At very high pressures, repulsive forces dominate and $Z > 1$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "T_fixed = 250 # K (above Tc for methane, supercritical)\n", + "P_range = np.logspace(4, 7.5, 100) # 10 kPa to ~30 MPa\n", + "\n", + "Z_pr_arr = []\n", + "Z_rks_arr = []\n", + "\n", + "for P in P_range:\n", + " _, Z = eval_eos(pr, T_fixed, P)\n", + " Z_pr_arr.append(Z)\n", + " _, Z = eval_eos(rks, T_fixed, P)\n", + " Z_rks_arr.append(Z)\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "ax.semilogx(P_range / 1e6, Z_pr_arr, label=\"Peng-Robinson\")\n", + "ax.semilogx(P_range / 1e6, Z_rks_arr, \"--\", label=\"Soave-Redlich-Kwong\")\n", + "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"Ideal gas\")\n", + "ax.set_xlabel(\"Pressure [MPa]\")\n", + "ax.set_ylabel(\"Compressibility Factor Z\")\n", + "ax.set_title(f\"Methane at T = {T_fixed} K\")\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Both EoS give very similar results. The characteristic dip below $Z = 1$ at moderate pressures reflects attractive intermolecular forces, while the rise above $Z = 1$ at high pressures is due to repulsive (excluded volume) effects." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Z vs Reduced Temperature\n", + "\n", + "Vary the temperature at a fixed pressure to see how the compressibility factor changes with reduced temperature $T_r = T / T_c$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P_fixed = 5e6 # 5 MPa\n", + "T_range = np.linspace(200, 600, 100)\n", + "Tr_range = T_range / Tc_CH4\n", + "\n", + "Z_vs_T = []\n", + "for T in T_range:\n", + " _, Z = eval_eos(pr, T, P_fixed)\n", + " Z_vs_T.append(Z)\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "ax.plot(Tr_range, Z_vs_T)\n", + "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5)\n", + "ax.axvline(1.0, color=\"red\", linestyle=\":\", alpha=0.5, label=r\"$T_r = 1$ (critical)\")\n", + "ax.set_xlabel(r\"Reduced Temperature $T_r = T / T_c$\")\n", + "ax.set_ylabel(\"Compressibility Factor Z\")\n", + "ax.set_title(f\"Peng-Robinson: Methane at P = {P_fixed/1e6:.0f} MPa\")\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Methane-Ethane Mixture\n", + "\n", + "The EoS blocks support mixtures through standard van der Waals one-fluid mixing rules. Supply arrays of critical properties and mole fractions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Critical properties: methane, ethane\n", + "Tc = [190.6, 305.3] # K\n", + "Pc = [4.6e6, 4.872e6] # Pa\n", + "omega = [0.011, 0.099] # acentric factors\n", + "\n", + "# Compare Z at different compositions\n", + "x1_range = np.linspace(0, 1, 20) # methane mole fraction\n", + "T_mix, P_mix = 300, 3e6 # 300 K, 3 MPa\n", + "\n", + "Z_mix = []\n", + "for x1 in x1_range:\n", + " pr_mix = PengRobinson(\n", + " Tc=Tc, Pc=Pc, omega=omega,\n", + " x=[x1, 1 - x1],\n", + " )\n", + " _, Z = eval_eos(pr_mix, T_mix, P_mix)\n", + " Z_mix.append(Z)\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "ax.plot(x1_range, Z_mix, \"o-\")\n", + "ax.set_xlabel(r\"$x_{\\mathrm{CH_4}}$\")\n", + "ax.set_ylabel(\"Compressibility Factor Z\")\n", + "ax.set_title(f\"PR: Methane-Ethane at T = {T_mix} K, P = {P_mix/1e6:.0f} MPa\")\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The mixture compressibility factor varies smoothly with composition. Pure ethane (right side) has a lower $Z$ at these conditions because it is closer to its critical point ($T_c = 305.3$ K) and therefore deviates more from ideal gas behavior." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/source/examples/vapor_pressure_curves.ipynb b/docs/source/examples/vapor_pressure_curves.ipynb new file mode 100644 index 0000000..c05aca8 --- /dev/null +++ b/docs/source/examples/vapor_pressure_curves.ipynb @@ -0,0 +1,195 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Vapor Pressure Curves\n", + "\n", + "Computing and comparing vapor pressure curves for water using different IK-CAPE correlation models: Antoine, Kirchhoff, and Wagner equations." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Import the correlation blocks and plotting libraries." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim_chem.thermodynamics import Antoine, Kirchhoff, Wagner" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Antoine Equation\n", + "\n", + "The Antoine equation is the most common three-parameter vapor pressure correlation:\n", + "\n", + "$$\\ln P_{\\text{sat}} = a_0 - \\frac{a_1}{T + a_2}$$\n", + "\n", + "We use well-known coefficients for water (natural log, T in Kelvin, P in Pa)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Antoine coefficients for water (NIST, natural log form)\n", + "antoine = Antoine(a0=23.2256, a1=3835.18, a2=-45.343)\n", + "\n", + "# Verify at the normal boiling point (373.15 K = 100 °C)\n", + "antoine.inputs[0] = 373.15\n", + "antoine.update(None)\n", + "P_sat = antoine.outputs[0]\n", + "print(f\"P_sat at 100 °C: {P_sat:.0f} Pa (expected ~101325 Pa)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Kirchhoff Equation\n", + "\n", + "The Kirchhoff equation is another three-parameter form:\n", + "\n", + "$$\\ln P_{\\text{sat}} = a_0 - \\frac{a_1}{T} - a_2 \\ln(T)$$\n", + "\n", + "It provides a slightly different functional form that can be more accurate over wider temperature ranges." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "# Kirchhoff coefficients for water (fitted to Antoine reference points)\nkirchhoff = Kirchhoff(a0=50.4338, a1=6350.4068, a2=3.6963)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Wagner Equation\n", + "\n", + "The Wagner equation uses reduced temperature and is considered one of the most accurate vapor pressure correlations:\n", + "\n", + "$$\\ln P_{\\text{sat}} = \\ln P_c + \\frac{1}{T_r}\\left(a_2 \\tau + a_3 \\tau^{1.5} + a_4 \\tau^{3} + a_5 \\tau^{6}\\right)$$\n", + "\n", + "where $T_r = T / T_c$ and $\\tau = 1 - T_r$. It requires the critical temperature and pressure as parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Wagner coefficients for water\n", + "# Tc = 647.096 K, Pc = 22064000 Pa\n", + "wagner = Wagner(\n", + " Tc=647.096, Pc=22064000,\n", + " a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Comparison\n", + "\n", + "Evaluate all three correlations over a temperature range and plot the vapor pressure curves. Each block follows the standard PathSim interface: set the input, call `update`, read the output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Temperature range: 300 K to 500 K\n", + "T_range = np.linspace(300, 500, 200)\n", + "\n", + "def eval_over_range(block, T_range):\n", + " \"\"\"Evaluate a correlation block over an array of temperatures.\"\"\"\n", + " results = []\n", + " for T in T_range:\n", + " block.inputs[0] = T\n", + " block.update(None)\n", + " results.append(block.outputs[0])\n", + " return np.array(results)\n", + "\n", + "P_antoine = eval_over_range(antoine, T_range)\n", + "P_kirchhoff = eval_over_range(kirchhoff, T_range)\n", + "P_wagner = eval_over_range(wagner, T_range)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + "\n", + "# Linear scale\n", + "ax1.plot(T_range - 273.15, P_antoine / 1e3, label=\"Antoine\")\n", + "ax1.plot(T_range - 273.15, P_kirchhoff / 1e3, \"--\", label=\"Kirchhoff\")\n", + "ax1.plot(T_range - 273.15, P_wagner / 1e3, \":\", label=\"Wagner\")\n", + "ax1.axhline(101.325, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"1 atm\")\n", + "ax1.set_xlabel(\"Temperature [°C]\")\n", + "ax1.set_ylabel(\"Vapor Pressure [kPa]\")\n", + "ax1.set_title(\"Vapor Pressure of Water\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Clausius-Clapeyron plot (log P vs 1/T)\n", + "ax2.semilogy(1000 / T_range, P_antoine, label=\"Antoine\")\n", + "ax2.semilogy(1000 / T_range, P_kirchhoff, \"--\", label=\"Kirchhoff\")\n", + "ax2.semilogy(1000 / T_range, P_wagner, \":\", label=\"Wagner\")\n", + "ax2.set_xlabel(\"1000 / T [1/K]\")\n", + "ax2.set_ylabel(\"Vapor Pressure [Pa]\")\n", + "ax2.set_title(\"Clausius-Clapeyron Plot\")\n", + "ax2.legend()\n", + "ax2.grid(True, alpha=0.3)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All three correlations agree well in their common validity range. The Wagner equation is typically preferred for high-accuracy work since it is constrained to reach the critical point." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/docs/source/examples/vle_calculation.ipynb b/docs/source/examples/vle_calculation.ipynb new file mode 100644 index 0000000..b587415 --- /dev/null +++ b/docs/source/examples/vle_calculation.ipynb @@ -0,0 +1,238 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# VLE with Activity Coefficients\n", + "\n", + "Computing vapor-liquid equilibrium (VLE) for the ethanol-water system using the NRTL activity coefficient model and Antoine vapor pressure correlation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Import the NRTL activity coefficient model and Antoine vapor pressure correlation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim_chem.thermodynamics import NRTL, Antoine" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pure Component Vapor Pressures\n", + "\n", + "First we set up Antoine correlations for both components. The modified Raoult's law gives:\n", + "\n", + "$$y_i P = x_i \\gamma_i P_i^{\\text{sat}}(T)$$\n", + "\n", + "At a fixed pressure, this determines the bubble-point temperature and vapor composition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Antoine coefficients (ln form, T in K, P in Pa)\n", + "# Ethanol\n", + "antoine_ethanol = Antoine(a0=23.5807, a1=3673.81, a2=-46.681)\n", + "# Water\n", + "antoine_water = Antoine(a0=23.2256, a1=3835.18, a2=-45.343)\n", + "\n", + "def psat(block, T):\n", + " \"\"\"Evaluate an Antoine block at temperature T.\"\"\"\n", + " block.inputs[0] = T\n", + " block.update(None)\n", + " return block.outputs[0]\n", + "\n", + "# Check normal boiling points\n", + "print(f\"Ethanol P_sat at 78.37 °C: {psat(antoine_ethanol, 351.52):.0f} Pa\")\n", + "print(f\"Water P_sat at 100 °C: {psat(antoine_water, 373.15):.0f} Pa\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## NRTL Activity Coefficients\n", + "\n", + "The NRTL model computes activity coefficients $\\gamma_1, \\gamma_2$ as a function of temperature and (fixed) composition. These quantify the deviation from ideal mixing." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_gammas(x1, T):\n", + " \"\"\"Compute NRTL activity coefficients for ethanol(1)-water(2).\"\"\"\n", + " nrtl = NRTL(\n", + " x=[x1, 1.0 - x1],\n", + " a=[[0, -0.801], [3.458, 0]],\n", + " c=[[0, 0.3], [0.3, 0]],\n", + " )\n", + " nrtl.inputs[0] = T\n", + " nrtl.update(None)\n", + " return nrtl.outputs[0], nrtl.outputs[1]\n", + "\n", + "# Activity coefficients at 50/50 mixture, 350 K\n", + "g1, g2 = get_gammas(0.5, 350)\n", + "print(f\"gamma_ethanol = {g1:.4f}\")\n", + "print(f\"gamma_water = {g2:.4f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Activity Coefficient vs Composition\n", + "\n", + "Plot how the activity coefficients vary with liquid composition at a fixed temperature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x1_range = np.linspace(0.01, 0.99, 50)\n", + "T_fixed = 351.0 # K, near ethanol boiling point\n", + "\n", + "gammas_1 = []\n", + "gammas_2 = []\n", + "for x1 in x1_range:\n", + " g1, g2 = get_gammas(x1, T_fixed)\n", + " gammas_1.append(g1)\n", + " gammas_2.append(g2)\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "ax.plot(x1_range, gammas_1, label=r\"$\\gamma_{\\mathrm{ethanol}}$\")\n", + "ax.plot(x1_range, gammas_2, label=r\"$\\gamma_{\\mathrm{water}}$\")\n", + "ax.axhline(1.0, color=\"gray\", linestyle=\"--\", alpha=0.5)\n", + "ax.set_xlabel(r\"$x_{\\mathrm{ethanol}}$\")\n", + "ax.set_ylabel(r\"$\\gamma_i$\")\n", + "ax.set_title(f\"NRTL Activity Coefficients at T = {T_fixed} K\")\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Both activity coefficients approach 1 as the mixture approaches pure component, as expected. The large values at infinite dilution (small $x_i$) reflect the strong non-ideal interactions in this system." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## T-x-y Diagram\n", + "\n", + "Construct the bubble-point T-x-y diagram at atmospheric pressure. For each liquid composition $x_1$, we solve for the temperature $T$ where:\n", + "\n", + "$$\\sum_i x_i \\gamma_i(T) P_i^{\\text{sat}}(T) = P$$\n", + "\n", + "and the vapor composition follows from:\n", + "\n", + "$$y_i = \\frac{x_i \\gamma_i P_i^{\\text{sat}}}{P}$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import brentq\n", + "\n", + "P_total = 101325 # 1 atm in Pa\n", + "\n", + "def bubble_pressure(x1, T):\n", + " \"\"\"Compute sum(x_i * gamma_i * P_sat_i) - P for bubble point.\"\"\"\n", + " x2 = 1.0 - x1\n", + " g1, g2 = get_gammas(x1, T)\n", + " P1 = psat(antoine_ethanol, T)\n", + " P2 = psat(antoine_water, T)\n", + " return x1 * g1 * P1 + x2 * g2 * P2 - P_total\n", + "\n", + "x1_points = np.linspace(0.01, 0.99, 40)\n", + "T_bubble = []\n", + "y1_values = []\n", + "\n", + "for x1 in x1_points:\n", + " # Solve for bubble-point T between 340 K and 380 K\n", + " T_bp = brentq(lambda T: bubble_pressure(x1, T), 340, 380)\n", + " T_bubble.append(T_bp)\n", + " \n", + " # Vapor composition\n", + " g1, g2 = get_gammas(x1, T_bp)\n", + " P1 = psat(antoine_ethanol, T_bp)\n", + " y1 = x1 * g1 * P1 / P_total\n", + " y1_values.append(y1)\n", + "\n", + "T_bubble = np.array(T_bubble)\n", + "y1_values = np.array(y1_values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "ax.plot(x1_points, T_bubble - 273.15, label=\"Bubble line (liquid)\")\n", + "ax.plot(y1_values, T_bubble - 273.15, label=\"Dew line (vapor)\")\n", + "ax.plot([0, 1], [0, 1], \"k--\", alpha=0.2, label=\"y = x\")\n", + "ax.set_xlabel(r\"$x_{\\mathrm{ethanol}}, \\ y_{\\mathrm{ethanol}}$\")\n", + "ax.set_ylabel(\"Temperature [°C]\")\n", + "ax.set_title(\"Ethanol-Water T-x-y Diagram at 1 atm\")\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The diagram shows the characteristic minimum-boiling azeotrope of the ethanol-water system near $x_{\\mathrm{ethanol}} \\approx 0.9$, which is a well-known feature of this system that prevents simple distillation from producing pure ethanol." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 784a68118d9ae1f289110879116f49805f852fb0 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 15:16:07 +0100 Subject: [PATCH 7/9] Add CI workflow for running tests on push and PR --- .github/workflows/tests.yml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .github/workflows/tests.yml diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..c81f204 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,18 @@ +name: Run tests +on: [push, pull_request] +jobs: + test: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.x' + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e .[test] + - name: Test with pytest + run: | + pytest tests/ -v From b9194d52b2c1ef1054c9050a3e4819f58c87ac22 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 15:18:15 +0100 Subject: [PATCH 8/9] Fix splitter test referencing removed _port_map_out attribute --- tests/tritium/test_splitter.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/tritium/test_splitter.py b/tests/tritium/test_splitter.py index 71f3e97..acc2013 100644 --- a/tests/tritium/test_splitter.py +++ b/tests/tritium/test_splitter.py @@ -35,8 +35,6 @@ def test_init(self): S = Splitter([0.4, 0.5, 0.1]) self.assertEqual(sum(S.fractions - np.array([0.4, 0.5, 0.1])), 0) - #test the automatic port maps - self.assertEqual(S._port_map_out, {"out 0.4":0, "out 0.5":1, "out 0.1":2}) def test_update(self): From 0c3f907f6fbbc7069db6cea2dc21835058265cf8 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Fri, 27 Feb 2026 15:25:14 +0100 Subject: [PATCH 9/9] Rewrite example notebooks to use Simulation and Connection --- docs/source/examples/equation_of_state.ipynb | 187 ++++++++-------- .../examples/vapor_pressure_curves.ipynb | 136 ++++-------- docs/source/examples/vle_calculation.ipynb | 208 +++++++----------- 3 files changed, 209 insertions(+), 322 deletions(-) diff --git a/docs/source/examples/equation_of_state.ipynb b/docs/source/examples/equation_of_state.ipynb index b8a9fd5..e04226d 100644 --- a/docs/source/examples/equation_of_state.ipynb +++ b/docs/source/examples/equation_of_state.ipynb @@ -6,16 +6,16 @@ "source": [ "# Cubic Equations of State\n", "\n", - "Computing compressibility factors and molar volumes using the Peng-Robinson and Soave-Redlich-Kwong cubic equations of state for pure methane and a methane-ethane mixture." + "Simulating the compressibility factor of methane across a pressure sweep using the Peng-Robinson and Soave-Redlich-Kwong equations of state wired into a PathSim simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Setup\n", + "Cubic equations of state compute the compressibility factor $Z = Pv/(RT)$ by solving a cubic polynomial at each $(T, P)$ condition. The EoS blocks take two inputs (temperature and pressure) and produce two outputs (molar volume and compressibility factor).\n", "\n", - "Import the equation of state blocks." + "At low pressure $Z \\to 1$ (ideal gas). At moderate pressure attractive forces cause $Z < 1$, and at high pressure excluded-volume effects push $Z > 1$." ] }, { @@ -24,9 +24,11 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Constant, Scope\n", + "\n", "from pathsim_chem.thermodynamics import PengRobinson, RedlichKwongSoave" ] }, @@ -34,13 +36,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Pure Methane\n", - "\n", - "The Peng-Robinson EoS solves a cubic equation in the compressibility factor $Z = Pv/(RT)$:\n", + "## System Definition\n", "\n", - "$$Z^3 - (1-B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0$$\n", - "\n", - "where $A$ and $B$ are dimensionless parameters computed from the critical properties. Each block takes $(T, P)$ as inputs and returns $(v, Z)$." + "Create Peng-Robinson and SRK blocks for pure methane. A `Constant` block supplies a fixed temperature while a `Source` sweeps pressure from 0.1 MPa to 30 MPa." ] }, { @@ -50,38 +48,29 @@ "outputs": [], "source": [ "# Critical properties of methane\n", - "Tc_CH4 = 190.6 # K\n", - "Pc_CH4 = 4.6e6 # Pa\n", - "omega_CH4 = 0.011 # acentric factor\n", - "\n", - "# Create EoS blocks\n", - "pr = PengRobinson(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", - "rks = RedlichKwongSoave(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", - "\n", - "def eval_eos(block, T, P):\n", - " \"\"\"Evaluate an EoS block and return (v, Z).\"\"\"\n", - " block.inputs[0] = T\n", - " block.inputs[1] = P\n", - " block.update(None)\n", - " return block.outputs[0], block.outputs[1]\n", - "\n", - "# Evaluate at 300 K, 1 atm\n", - "v_pr, Z_pr = eval_eos(pr, 300, 101325)\n", - "v_rks, Z_rks = eval_eos(rks, 300, 101325)\n", - "\n", - "print(f\"Methane at 300 K, 1 atm:\")\n", - "print(f\" PR: Z = {Z_pr:.6f}, v = {v_pr:.6e} m³/mol\")\n", - "print(f\" RKS: Z = {Z_rks:.6f}, v = {v_rks:.6e} m³/mol\")\n", - "print(f\" Ideal gas: Z = 1.0, v = {8.314462 * 300 / 101325:.6e} m³/mol\")" + "Tc, Pc, omega = 190.6, 4.6e6, 0.011\n", + "\n", + "# EoS blocks\n", + "pr = PengRobinson(Tc=Tc, Pc=Pc, omega=omega)\n", + "rks = RedlichKwongSoave(Tc=Tc, Pc=Pc, omega=omega)\n", + "\n", + "# Fixed temperature, logarithmic pressure sweep\n", + "import numpy as np\n", + "T_const = Constant(250) # 250 K (above Tc, supercritical)\n", + "P_src = Source(func=lambda t: 10**(4 + t * 0.035)) # 10 kPa to ~30 MPa over 100s\n", + "\n", + "# Scopes: record Z from both EoS (output port 1)\n", + "scp_pr = Scope(labels=[\"Z_PR\"])\n", + "scp_rks = Scope(labels=[\"Z_RKS\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Compressibility Factor vs Pressure\n", + "## Wiring\n", "\n", - "At low pressure, $Z \\to 1$ (ideal gas). As pressure increases, intermolecular forces cause $Z$ to deviate. At very high pressures, repulsive forces dominate and $Z > 1$." + "Both EoS blocks receive the same $(T, P)$ inputs. We record the compressibility factor (output port 1) from each." ] }, { @@ -90,25 +79,40 @@ "metadata": {}, "outputs": [], "source": [ - "T_fixed = 250 # K (above Tc for methane, supercritical)\n", - "P_range = np.logspace(4, 7.5, 100) # 10 kPa to ~30 MPa\n", - "\n", - "Z_pr_arr = []\n", - "Z_rks_arr = []\n", - "\n", - "for P in P_range:\n", - " _, Z = eval_eos(pr, T_fixed, P)\n", - " Z_pr_arr.append(Z)\n", - " _, Z = eval_eos(rks, T_fixed, P)\n", - " Z_rks_arr.append(Z)\n", + "sim = Simulation(\n", + " blocks=[T_const, P_src, pr, rks, scp_pr, scp_rks],\n", + " connections=[\n", + " # Temperature -> both EoS (input port 0)\n", + " Connection(T_const, pr, rks),\n", + " # Pressure -> both EoS (input port 1)\n", + " Connection(P_src, pr[1], rks[1]),\n", + " # Z output (port 1) -> scopes\n", + " Connection(pr[1], scp_pr),\n", + " Connection(rks[1], scp_rks),\n", + " ],\n", + " dt=1.0,\n", + ")\n", + "\n", + "sim.run(100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time, Z_pr = scp_pr.read()\n", + "_, Z_rks = scp_rks.read()\n", + "P_vals = 10**(4 + time * 0.035) / 1e6 # MPa\n", "\n", "fig, ax = plt.subplots(figsize=(7, 5))\n", - "ax.semilogx(P_range / 1e6, Z_pr_arr, label=\"Peng-Robinson\")\n", - "ax.semilogx(P_range / 1e6, Z_rks_arr, \"--\", label=\"Soave-Redlich-Kwong\")\n", + "ax.semilogx(P_vals, Z_pr[0], label=\"Peng-Robinson\")\n", + "ax.semilogx(P_vals, Z_rks[0], \"--\", label=\"Soave-Redlich-Kwong\")\n", "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"Ideal gas\")\n", "ax.set_xlabel(\"Pressure [MPa]\")\n", "ax.set_ylabel(\"Compressibility Factor Z\")\n", - "ax.set_title(f\"Methane at T = {T_fixed} K\")\n", + "ax.set_title(\"Methane at T = 250 K\")\n", "ax.legend()\n", "ax.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", @@ -126,9 +130,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Z vs Reduced Temperature\n", + "## Mixture\n", "\n", - "Vary the temperature at a fixed pressure to see how the compressibility factor changes with reduced temperature $T_r = T / T_c$." + "The EoS blocks also support mixtures through van der Waals one-fluid mixing rules. Here we set up a methane-ethane mixture and sweep pressure." ] }, { @@ -137,35 +141,28 @@ "metadata": {}, "outputs": [], "source": [ - "P_fixed = 5e6 # 5 MPa\n", - "T_range = np.linspace(200, 600, 100)\n", - "Tr_range = T_range / Tc_CH4\n", - "\n", - "Z_vs_T = []\n", - "for T in T_range:\n", - " _, Z = eval_eos(pr, T, P_fixed)\n", - " Z_vs_T.append(Z)\n", - "\n", - "fig, ax = plt.subplots(figsize=(7, 5))\n", - "ax.plot(Tr_range, Z_vs_T)\n", - "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5)\n", - "ax.axvline(1.0, color=\"red\", linestyle=\":\", alpha=0.5, label=r\"$T_r = 1$ (critical)\")\n", - "ax.set_xlabel(r\"Reduced Temperature $T_r = T / T_c$\")\n", - "ax.set_ylabel(\"Compressibility Factor Z\")\n", - "ax.set_title(f\"Peng-Robinson: Methane at P = {P_fixed/1e6:.0f} MPa\")\n", - "ax.legend()\n", - "ax.grid(True, alpha=0.3)\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Methane-Ethane Mixture\n", - "\n", - "The EoS blocks support mixtures through standard van der Waals one-fluid mixing rules. Supply arrays of critical properties and mole fractions." + "pr_mix = PengRobinson(\n", + " Tc=[190.6, 305.3],\n", + " Pc=[4.6e6, 4.872e6],\n", + " omega=[0.011, 0.099],\n", + " x=[0.7, 0.3],\n", + ")\n", + "\n", + "T_const2 = Constant(300)\n", + "P_src2 = Source(func=lambda t: 10**(4 + t * 0.035))\n", + "scp_mix = Scope(labels=[\"Z_mixture\"])\n", + "\n", + "sim_mix = Simulation(\n", + " blocks=[T_const2, P_src2, pr_mix, scp_mix],\n", + " connections=[\n", + " Connection(T_const2, pr_mix),\n", + " Connection(P_src2, pr_mix[1]),\n", + " Connection(pr_mix[1], scp_mix),\n", + " ],\n", + " dt=1.0,\n", + ")\n", + "\n", + "sim_mix.run(100)" ] }, { @@ -174,29 +171,17 @@ "metadata": {}, "outputs": [], "source": [ - "# Critical properties: methane, ethane\n", - "Tc = [190.6, 305.3] # K\n", - "Pc = [4.6e6, 4.872e6] # Pa\n", - "omega = [0.011, 0.099] # acentric factors\n", - "\n", - "# Compare Z at different compositions\n", - "x1_range = np.linspace(0, 1, 20) # methane mole fraction\n", - "T_mix, P_mix = 300, 3e6 # 300 K, 3 MPa\n", - "\n", - "Z_mix = []\n", - "for x1 in x1_range:\n", - " pr_mix = PengRobinson(\n", - " Tc=Tc, Pc=Pc, omega=omega,\n", - " x=[x1, 1 - x1],\n", - " )\n", - " _, Z = eval_eos(pr_mix, T_mix, P_mix)\n", - " Z_mix.append(Z)\n", + "time_m, Z_mix = scp_mix.read()\n", + "P_mix = 10**(4 + time_m * 0.035) / 1e6\n", "\n", "fig, ax = plt.subplots(figsize=(7, 5))\n", - "ax.plot(x1_range, Z_mix, \"o-\")\n", - "ax.set_xlabel(r\"$x_{\\mathrm{CH_4}}$\")\n", + "ax.semilogx(P_vals, Z_pr[0], label=\"Pure CH₄ (250 K)\")\n", + "ax.semilogx(P_mix, Z_mix[0], \"--\", label=\"70/30 CH₄-C₂H₆ (300 K)\")\n", + "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5)\n", + "ax.set_xlabel(\"Pressure [MPa]\")\n", "ax.set_ylabel(\"Compressibility Factor Z\")\n", - "ax.set_title(f\"PR: Methane-Ethane at T = {T_mix} K, P = {P_mix/1e6:.0f} MPa\")\n", + "ax.set_title(\"Peng-Robinson: Pure vs Mixture\")\n", + "ax.legend()\n", "ax.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "plt.show()" @@ -206,7 +191,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The mixture compressibility factor varies smoothly with composition. Pure ethane (right side) has a lower $Z$ at these conditions because it is closer to its critical point ($T_c = 305.3$ K) and therefore deviates more from ideal gas behavior." + "The mixture shows a deeper dip because ethane ($T_c = 305.3$ K) is near its critical temperature at 300 K, leading to stronger non-ideal behavior." ] } ], diff --git a/docs/source/examples/vapor_pressure_curves.ipynb b/docs/source/examples/vapor_pressure_curves.ipynb index c05aca8..0191f36 100644 --- a/docs/source/examples/vapor_pressure_curves.ipynb +++ b/docs/source/examples/vapor_pressure_curves.ipynb @@ -6,16 +6,20 @@ "source": [ "# Vapor Pressure Curves\n", "\n", - "Computing and comparing vapor pressure curves for water using different IK-CAPE correlation models: Antoine, Kirchhoff, and Wagner equations." + "Comparing vapor pressure correlations for water by sweeping temperature through Antoine, Kirchhoff, and Wagner blocks wired into a PathSim simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Setup\n", + "The IK-CAPE standard defines several temperature-dependent correlations for vapor pressure. Here we compare three of them:\n", "\n", - "Import the correlation blocks and plotting libraries." + "- **Antoine**: $\\ln P = a_0 - a_1 / (T + a_2)$\n", + "- **Kirchhoff**: $\\ln P = a_0 - a_1/T - a_2 \\ln T$\n", + "- **Wagner**: $\\ln P = \\ln P_c + (1/T_r)(a_2\\tau + a_3\\tau^{1.5} + a_4\\tau^3 + a_5\\tau^6)$\n", + "\n", + "Each correlation is a `Function` block with a single temperature input and a single pressure output." ] }, { @@ -24,9 +28,11 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Scope\n", + "\n", "from pathsim_chem.thermodynamics import Antoine, Kirchhoff, Wagner" ] }, @@ -34,13 +40,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Antoine Equation\n", - "\n", - "The Antoine equation is the most common three-parameter vapor pressure correlation:\n", - "\n", - "$$\\ln P_{\\text{sat}} = a_0 - \\frac{a_1}{T + a_2}$$\n", - "\n", - "We use well-known coefficients for water (natural log, T in Kelvin, P in Pa)." + "Define the correlation blocks with literature coefficients for water. All use natural logarithm, temperature in Kelvin, and return pressure in Pascals." ] }, { @@ -49,70 +49,21 @@ "metadata": {}, "outputs": [], "source": [ - "# Antoine coefficients for water (NIST, natural log form)\n", + "# Antoine coefficients for water (NIST)\n", "antoine = Antoine(a0=23.2256, a1=3835.18, a2=-45.343)\n", "\n", - "# Verify at the normal boiling point (373.15 K = 100 °C)\n", - "antoine.inputs[0] = 373.15\n", - "antoine.update(None)\n", - "P_sat = antoine.outputs[0]\n", - "print(f\"P_sat at 100 °C: {P_sat:.0f} Pa (expected ~101325 Pa)\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Kirchhoff Equation\n", - "\n", - "The Kirchhoff equation is another three-parameter form:\n", - "\n", - "$$\\ln P_{\\text{sat}} = a_0 - \\frac{a_1}{T} - a_2 \\ln(T)$$\n", - "\n", - "It provides a slightly different functional form that can be more accurate over wider temperature ranges." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": "# Kirchhoff coefficients for water (fitted to Antoine reference points)\nkirchhoff = Kirchhoff(a0=50.4338, a1=6350.4068, a2=3.6963)" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Wagner Equation\n", - "\n", - "The Wagner equation uses reduced temperature and is considered one of the most accurate vapor pressure correlations:\n", - "\n", - "$$\\ln P_{\\text{sat}} = \\ln P_c + \\frac{1}{T_r}\\left(a_2 \\tau + a_3 \\tau^{1.5} + a_4 \\tau^{3} + a_5 \\tau^{6}\\right)$$\n", + "# Kirchhoff coefficients for water\n", + "kirchhoff = Kirchhoff(a0=50.4338, a1=6350.4068, a2=3.6963)\n", "\n", - "where $T_r = T / T_c$ and $\\tau = 1 - T_r$. It requires the critical temperature and pressure as parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Wagner coefficients for water\n", - "# Tc = 647.096 K, Pc = 22064000 Pa\n", - "wagner = Wagner(\n", - " Tc=647.096, Pc=22064000,\n", - " a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303,\n", - ")" + "# Wagner coefficients for water (Tc=647.096 K, Pc=22.064 MPa)\n", + "wagner = Wagner(Tc=647.096, Pc=22064000, a2=-7.76451, a3=1.45838, a4=-2.77580, a5=-1.23303)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Comparison\n", - "\n", - "Evaluate all three correlations over a temperature range and plot the vapor pressure curves. Each block follows the standard PathSim interface: set the input, call `update`, read the output." + "Set up a simulation that sweeps temperature from 300 K to 500 K. A `Source` block generates a linear temperature ramp, which is connected to all three correlations in parallel. A `Scope` records the three pressure outputs." ] }, { @@ -121,21 +72,24 @@ "metadata": {}, "outputs": [], "source": [ - "# Temperature range: 300 K to 500 K\n", - "T_range = np.linspace(300, 500, 200)\n", - "\n", - "def eval_over_range(block, T_range):\n", - " \"\"\"Evaluate a correlation block over an array of temperatures.\"\"\"\n", - " results = []\n", - " for T in T_range:\n", - " block.inputs[0] = T\n", - " block.update(None)\n", - " results.append(block.outputs[0])\n", - " return np.array(results)\n", - "\n", - "P_antoine = eval_over_range(antoine, T_range)\n", - "P_kirchhoff = eval_over_range(kirchhoff, T_range)\n", - "P_wagner = eval_over_range(wagner, T_range)" + "# Temperature ramp: 300 K to 500 K over 200 seconds\n", + "T_src = Source(func=lambda t: 300 + t)\n", + "\n", + "# Scope to record all three vapor pressures\n", + "scp = Scope(labels=[\"Antoine\", \"Kirchhoff\", \"Wagner\"])\n", + "\n", + "sim = Simulation(\n", + " blocks=[T_src, antoine, kirchhoff, wagner, scp],\n", + " connections=[\n", + " Connection(T_src, antoine, kirchhoff, wagner),\n", + " Connection(antoine, scp),\n", + " Connection(kirchhoff, scp[1]),\n", + " Connection(wagner, scp[2]),\n", + " ],\n", + " dt=1.0,\n", + ")\n", + "\n", + "sim.run(200)" ] }, { @@ -144,12 +98,18 @@ "metadata": {}, "outputs": [], "source": [ + "time, signals = scp.read()\n", + "\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", "\n", - "# Linear scale\n", - "ax1.plot(T_range - 273.15, P_antoine / 1e3, label=\"Antoine\")\n", - "ax1.plot(T_range - 273.15, P_kirchhoff / 1e3, \"--\", label=\"Kirchhoff\")\n", - "ax1.plot(T_range - 273.15, P_wagner / 1e3, \":\", label=\"Wagner\")\n", + "T_celsius = time + 300 - 273.15\n", + "labels = [\"Antoine\", \"Kirchhoff\", \"Wagner\"]\n", + "styles = [\"-\", \"--\", \":\"]\n", + "\n", + "for i, (label, style) in enumerate(zip(labels, styles)):\n", + " ax1.plot(T_celsius, signals[i] / 1e3, style, label=label)\n", + " ax2.semilogy(1000 / (time + 300), signals[i], style, label=label)\n", + "\n", "ax1.axhline(101.325, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"1 atm\")\n", "ax1.set_xlabel(\"Temperature [°C]\")\n", "ax1.set_ylabel(\"Vapor Pressure [kPa]\")\n", @@ -157,10 +117,6 @@ "ax1.legend()\n", "ax1.grid(True, alpha=0.3)\n", "\n", - "# Clausius-Clapeyron plot (log P vs 1/T)\n", - "ax2.semilogy(1000 / T_range, P_antoine, label=\"Antoine\")\n", - "ax2.semilogy(1000 / T_range, P_kirchhoff, \"--\", label=\"Kirchhoff\")\n", - "ax2.semilogy(1000 / T_range, P_wagner, \":\", label=\"Wagner\")\n", "ax2.set_xlabel(\"1000 / T [1/K]\")\n", "ax2.set_ylabel(\"Vapor Pressure [Pa]\")\n", "ax2.set_title(\"Clausius-Clapeyron Plot\")\n", @@ -192,4 +148,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} diff --git a/docs/source/examples/vle_calculation.ipynb b/docs/source/examples/vle_calculation.ipynb index b587415..99bf0cd 100644 --- a/docs/source/examples/vle_calculation.ipynb +++ b/docs/source/examples/vle_calculation.ipynb @@ -4,18 +4,20 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# VLE with Activity Coefficients\n", + "# Activity Coefficients\n", "\n", - "Computing vapor-liquid equilibrium (VLE) for the ethanol-water system using the NRTL activity coefficient model and Antoine vapor pressure correlation." + "Simulating how NRTL activity coefficients for the ethanol-water system evolve with temperature using PathSim blocks and connections." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Setup\n", + "The NRTL (Non-Random Two-Liquid) model computes liquid-phase activity coefficients $\\gamma_i$ that quantify deviations from ideal mixing:\n", "\n", - "Import the NRTL activity coefficient model and Antoine vapor pressure correlation." + "$$y_i P = x_i \\gamma_i P_i^{\\text{sat}}(T)$$\n", + "\n", + "Here we wire an NRTL block and two Antoine blocks into a simulation to see how activity coefficients and the resulting bubble pressures change with temperature." ] }, { @@ -24,9 +26,11 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Scope, Function\n", + "\n", "from pathsim_chem.thermodynamics import NRTL, Antoine" ] }, @@ -34,13 +38,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Pure Component Vapor Pressures\n", + "## System Definition\n", "\n", - "First we set up Antoine correlations for both components. The modified Raoult's law gives:\n", - "\n", - "$$y_i P = x_i \\gamma_i P_i^{\\text{sat}}(T)$$\n", - "\n", - "At a fixed pressure, this determines the bubble-point temperature and vapor composition." + "Set up the blocks: a temperature source, Antoine correlations for both components, and an NRTL activity coefficient model for a 40/60 ethanol-water mixture." ] }, { @@ -49,30 +49,29 @@ "metadata": {}, "outputs": [], "source": [ - "# Antoine coefficients (ln form, T in K, P in Pa)\n", - "# Ethanol\n", + "# Temperature ramp from 330 K to 380 K\n", + "T_src = Source(func=lambda t: 330 + t * 0.5)\n", + "\n", + "# Antoine vapor pressure correlations (ln form, T in K, P in Pa)\n", "antoine_ethanol = Antoine(a0=23.5807, a1=3673.81, a2=-46.681)\n", - "# Water\n", "antoine_water = Antoine(a0=23.2256, a1=3835.18, a2=-45.343)\n", "\n", - "def psat(block, T):\n", - " \"\"\"Evaluate an Antoine block at temperature T.\"\"\"\n", - " block.inputs[0] = T\n", - " block.update(None)\n", - " return block.outputs[0]\n", - "\n", - "# Check normal boiling points\n", - "print(f\"Ethanol P_sat at 78.37 °C: {psat(antoine_ethanol, 351.52):.0f} Pa\")\n", - "print(f\"Water P_sat at 100 °C: {psat(antoine_water, 373.15):.0f} Pa\")" + "# NRTL activity coefficients for ethanol(1)-water(2)\n", + "nrtl = NRTL(\n", + " x=[0.4, 0.6],\n", + " a=[[0, -0.801], [3.458, 0]],\n", + " b=[[0, 200], [-100, 0]],\n", + " c=[[0, 0.3], [0.3, 0]],\n", + ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## NRTL Activity Coefficients\n", + "A custom `Function` block computes the bubble pressure from the activity coefficients and pure-component vapor pressures:\n", "\n", - "The NRTL model computes activity coefficients $\\gamma_1, \\gamma_2$ as a function of temperature and (fixed) composition. These quantify the deviation from ideal mixing." + "$$P_{\\text{bubble}} = x_1 \\gamma_1 P_1^{\\text{sat}} + x_2 \\gamma_2 P_2^{\\text{sat}}$$" ] }, { @@ -81,81 +80,25 @@ "metadata": {}, "outputs": [], "source": [ - "def get_gammas(x1, T):\n", - " \"\"\"Compute NRTL activity coefficients for ethanol(1)-water(2).\"\"\"\n", - " nrtl = NRTL(\n", - " x=[x1, 1.0 - x1],\n", - " a=[[0, -0.801], [3.458, 0]],\n", - " c=[[0, 0.3], [0.3, 0]],\n", - " )\n", - " nrtl.inputs[0] = T\n", - " nrtl.update(None)\n", - " return nrtl.outputs[0], nrtl.outputs[1]\n", - "\n", - "# Activity coefficients at 50/50 mixture, 350 K\n", - "g1, g2 = get_gammas(0.5, 350)\n", - "print(f\"gamma_ethanol = {g1:.4f}\")\n", - "print(f\"gamma_water = {g2:.4f}\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Activity Coefficient vs Composition\n", + "# Compute bubble pressure from gamma and Psat values\n", + "x1, x2 = 0.4, 0.6\n", "\n", - "Plot how the activity coefficients vary with liquid composition at a fixed temperature." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x1_range = np.linspace(0.01, 0.99, 50)\n", - "T_fixed = 351.0 # K, near ethanol boiling point\n", - "\n", - "gammas_1 = []\n", - "gammas_2 = []\n", - "for x1 in x1_range:\n", - " g1, g2 = get_gammas(x1, T_fixed)\n", - " gammas_1.append(g1)\n", - " gammas_2.append(g2)\n", - "\n", - "fig, ax = plt.subplots(figsize=(7, 5))\n", - "ax.plot(x1_range, gammas_1, label=r\"$\\gamma_{\\mathrm{ethanol}}$\")\n", - "ax.plot(x1_range, gammas_2, label=r\"$\\gamma_{\\mathrm{water}}$\")\n", - "ax.axhline(1.0, color=\"gray\", linestyle=\"--\", alpha=0.5)\n", - "ax.set_xlabel(r\"$x_{\\mathrm{ethanol}}$\")\n", - "ax.set_ylabel(r\"$\\gamma_i$\")\n", - "ax.set_title(f\"NRTL Activity Coefficients at T = {T_fixed} K\")\n", - "ax.legend()\n", - "ax.grid(True, alpha=0.3)\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Both activity coefficients approach 1 as the mixture approaches pure component, as expected. The large values at infinite dilution (small $x_i$) reflect the strong non-ideal interactions in this system." + "bubble = Function(\n", + " func=lambda g1, g2, P1, P2: x1 * g1 * P1 + x2 * g2 * P2\n", + ")\n", + "\n", + "# Scopes to record activity coefficients and bubble pressure\n", + "scp_gamma = Scope(labels=[\"gamma_ethanol\", \"gamma_water\"])\n", + "scp_bubble = Scope(labels=[\"P_bubble\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## T-x-y Diagram\n", + "## Wiring and Simulation\n", "\n", - "Construct the bubble-point T-x-y diagram at atmospheric pressure. For each liquid composition $x_1$, we solve for the temperature $T$ where:\n", - "\n", - "$$\\sum_i x_i \\gamma_i(T) P_i^{\\text{sat}}(T) = P$$\n", - "\n", - "and the vapor composition follows from:\n", - "\n", - "$$y_i = \\frac{x_i \\gamma_i P_i^{\\text{sat}}}{P}$$" + "Connect the temperature source to NRTL and both Antoine blocks. The NRTL outputs and Antoine outputs feed into the bubble pressure calculation. Scopes record everything." ] }, { @@ -164,35 +107,24 @@ "metadata": {}, "outputs": [], "source": [ - "from scipy.optimize import brentq\n", - "\n", - "P_total = 101325 # 1 atm in Pa\n", - "\n", - "def bubble_pressure(x1, T):\n", - " \"\"\"Compute sum(x_i * gamma_i * P_sat_i) - P for bubble point.\"\"\"\n", - " x2 = 1.0 - x1\n", - " g1, g2 = get_gammas(x1, T)\n", - " P1 = psat(antoine_ethanol, T)\n", - " P2 = psat(antoine_water, T)\n", - " return x1 * g1 * P1 + x2 * g2 * P2 - P_total\n", - "\n", - "x1_points = np.linspace(0.01, 0.99, 40)\n", - "T_bubble = []\n", - "y1_values = []\n", - "\n", - "for x1 in x1_points:\n", - " # Solve for bubble-point T between 340 K and 380 K\n", - " T_bp = brentq(lambda T: bubble_pressure(x1, T), 340, 380)\n", - " T_bubble.append(T_bp)\n", - " \n", - " # Vapor composition\n", - " g1, g2 = get_gammas(x1, T_bp)\n", - " P1 = psat(antoine_ethanol, T_bp)\n", - " y1 = x1 * g1 * P1 / P_total\n", - " y1_values.append(y1)\n", - "\n", - "T_bubble = np.array(T_bubble)\n", - "y1_values = np.array(y1_values)" + "sim = Simulation(\n", + " blocks=[T_src, nrtl, antoine_ethanol, antoine_water, bubble, scp_gamma, scp_bubble],\n", + " connections=[\n", + " # Temperature drives all thermodynamic blocks\n", + " Connection(T_src, nrtl, antoine_ethanol, antoine_water),\n", + " # NRTL outputs -> scope and bubble pressure\n", + " Connection(nrtl, scp_gamma, bubble),\n", + " Connection(nrtl[1], scp_gamma[1], bubble[1]),\n", + " # Antoine outputs -> bubble pressure\n", + " Connection(antoine_ethanol, bubble[2]),\n", + " Connection(antoine_water, bubble[3]),\n", + " # Bubble pressure -> scope\n", + " Connection(bubble, scp_bubble),\n", + " ],\n", + " dt=0.5,\n", + ")\n", + "\n", + "sim.run(100)" ] }, { @@ -201,15 +133,29 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(figsize=(7, 5))\n", - "ax.plot(x1_points, T_bubble - 273.15, label=\"Bubble line (liquid)\")\n", - "ax.plot(y1_values, T_bubble - 273.15, label=\"Dew line (vapor)\")\n", - "ax.plot([0, 1], [0, 1], \"k--\", alpha=0.2, label=\"y = x\")\n", - "ax.set_xlabel(r\"$x_{\\mathrm{ethanol}}, \\ y_{\\mathrm{ethanol}}$\")\n", - "ax.set_ylabel(\"Temperature [°C]\")\n", - "ax.set_title(\"Ethanol-Water T-x-y Diagram at 1 atm\")\n", - "ax.legend()\n", - "ax.grid(True, alpha=0.3)\n", + "time_g, gamma = scp_gamma.read()\n", + "time_b, P_bub = scp_bubble.read()\n", + "T_range = 330 + time_g * 0.5\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + "\n", + "ax1.plot(T_range - 273.15, gamma[0], label=r\"$\\gamma_{\\mathrm{ethanol}}$\")\n", + "ax1.plot(T_range - 273.15, gamma[1], label=r\"$\\gamma_{\\mathrm{water}}$\")\n", + "ax1.axhline(1.0, color=\"gray\", linestyle=\"--\", alpha=0.5)\n", + "ax1.set_xlabel(\"Temperature [°C]\")\n", + "ax1.set_ylabel(r\"$\\gamma_i$\")\n", + "ax1.set_title(\"NRTL Activity Coefficients\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "ax2.plot(T_range - 273.15, P_bub[0] / 1e3)\n", + "ax2.axhline(101.325, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"1 atm\")\n", + "ax2.set_xlabel(\"Temperature [°C]\")\n", + "ax2.set_ylabel(\"Bubble Pressure [kPa]\")\n", + "ax2.set_title(\"Bubble Pressure (40/60 Ethanol-Water)\")\n", + "ax2.legend()\n", + "ax2.grid(True, alpha=0.3)\n", + "\n", "plt.tight_layout()\n", "plt.show()" ] @@ -218,7 +164,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The diagram shows the characteristic minimum-boiling azeotrope of the ethanol-water system near $x_{\\mathrm{ethanol}} \\approx 0.9$, which is a well-known feature of this system that prevents simple distillation from producing pure ethanol." + "Both activity coefficients exceed 1, confirming positive deviation from Raoult's law — a characteristic of the ethanol-water system. The bubble pressure rises steeply with temperature, crossing 1 atm near the expected boiling point of the mixture." ] } ],