diff --git a/docs/source/examples/cstr_reaction.ipynb b/docs/source/examples/cstr_reaction.ipynb new file mode 100644 index 0000000..918a9d6 --- /dev/null +++ b/docs/source/examples/cstr_reaction.ipynb @@ -0,0 +1,112 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Continuous Stirred-Tank Reactor\n", + "\n", + "Simulating the startup transient of an exothermic first-order reaction in a cooled CSTR, showing the dynamic interaction between concentration decay and temperature rise." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The CSTR model couples a material balance with an energy balance. For a first-order irreversible reaction $A \\to \\text{products}$ with Arrhenius kinetics:\n", + "\n", + "$$\\frac{dC_A}{dt} = \\frac{C_{A,\\text{in}} - C_A}{\\tau} - k(T)\\, C_A$$\n", + "\n", + "$$\\frac{dT}{dt} = \\frac{T_\\text{in} - T}{\\tau} + \\frac{(-\\Delta H_\\text{rxn})}{\\rho\\, C_p}\\, k(T)\\, C_A + \\frac{UA}{\\rho\\, C_p\\, V}\\,(T_c - T)$$\n", + "\n", + "where $k(T) = k_0 \\exp(-E_a / RT)$ and $\\tau = V/F$ is the residence time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Scope\n", + "\n", + "from pathsim_chem.process import CSTR" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Configure the reactor with parameters typical of an exothermic liquid-phase reaction. The reactor starts empty ($C_A = 0$) at ambient temperature and is fed with a concentrated stream." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "cstr = CSTR(\n V=1.0, # reactor volume [m³]\n F=0.05, # volumetric flow rate [m³/s] -> tau = 20 s\n k0=1e6, # pre-exponential factor [1/s]\n Ea=40000.0, # activation energy [J/mol]\n n=1.0, # first-order reaction\n dH_rxn=-40000.0, # exothermic [J/mol]\n rho=1000.0, # density [kg/m³]\n Cp=4184.0, # heat capacity [J/(kg·K)]\n UA=800.0, # cooling jacket [W/K]\n C_A0=0.0, # start empty\n T0=300.0, # initial temperature [K]\n)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "Feed a constant concentration of 1000 mol/m³ at 320 K, with the coolant held at 290 K. A `Scope` records both the outlet concentration and temperature." + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "# Constant feed and coolant conditions\nC_feed = Source(func=lambda t: 1000.0) # feed concentration [mol/m³]\nT_feed = Source(func=lambda t: 320.0) # feed temperature [K]\nT_cool = Source(func=lambda t: 290.0) # coolant temperature [K]\n\nscp = Scope(labels=[\"C_A [mol/m³]\", \"T [K]\"])\n\nsim = Simulation(\n blocks=[C_feed, T_feed, T_cool, cstr, scp],\n connections=[\n Connection(C_feed, cstr), # C_in -> port 0\n Connection(T_feed, cstr[1]), # T_in -> port 1\n Connection(T_cool, cstr[2]), # T_c -> port 2\n Connection(cstr, scp), # C_out -> scope port 0\n Connection(cstr[1], scp[1]), # T_out -> scope port 1\n ],\n dt=0.1,\n)\n\nsim.run(200)" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time, signals = scp.read()\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + "\n", + "ax1.plot(time, signals[0])\n", + "ax1.set_xlabel(\"Time [s]\")\n", + "ax1.set_ylabel(\"Concentration [mol/m³]\")\n", + "ax1.set_title(\"Outlet Concentration\")\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "ax2.plot(time, signals[1])\n", + "ax2.set_xlabel(\"Time [s]\")\n", + "ax2.set_ylabel(\"Temperature [K]\")\n", + "ax2.set_title(\"Reactor Temperature\")\n", + "ax2.grid(True, alpha=0.3)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The reactor starts cold and empty. As fresh feed enters, concentration rises initially but then the Arrhenius kinetics kick in — the exothermic reaction heats the fluid, which accelerates the rate, consuming more reactant. The cooling jacket prevents thermal runaway and the system settles to a steady state where reaction rate balances the feed." + ] + } + ], + "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/flash_distillation.ipynb b/docs/source/examples/flash_distillation.ipynb new file mode 100644 index 0000000..5e60cdd --- /dev/null +++ b/docs/source/examples/flash_distillation.ipynb @@ -0,0 +1,206 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Flash Drum and Distillation Column\n", + "\n", + "Simulating two fundamental separation processes: an isothermal binary flash drum and a multi-tray distillation column built from individual `DistillationTray` blocks wired in series." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Part 1: Isothermal Flash Drum\n", + "\n", + "A flash drum separates a liquid feed into vapor and liquid streams using vapor-liquid equilibrium (VLE). The drum uses Raoult's law with Antoine correlations to compute K-values:\n", + "\n", + "$$K_i = \\frac{P^\\text{sat}_i(T)}{P}$$\n", + "\n", + "The Rachford-Rice equation determines the vapor fraction $\\beta$, from which the vapor ($y_i$) and liquid ($x_i$) compositions follow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Scope\n", + "\n", + "from pathsim_chem.process import FlashDrum, DistillationTray" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "Configure a flash drum for a benzene-toluene mixture (default Antoine coefficients). Feed an equimolar mixture at 1 atm and sweep temperature from 340 K to 400 K. This range covers the bubble point (~365 K) and dew point (~380 K) of the mixture, so we can observe the transition from all-liquid to two-phase to all-vapor." + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "flash = FlashDrum(holdup=100.0) # default benzene/toluene Antoine params\n\n# Feed: 10 mol/s, equimolar (z1 = 0.5), 1 atm\nF_src = Source(func=lambda t: 10.0)\nz_src = Source(func=lambda t: 0.5)\nT_src = Source(func=lambda t: 340.0 + t * 0.5) # ramp 340 -> 400 K\nP_src = Source(func=lambda t: 101325.0)\n\nscp = Scope(labels=[\"V_rate\", \"L_rate\", \"y_1 (benzene)\", \"x_1 (benzene)\"])\n\nsim = Simulation(\n blocks=[F_src, z_src, T_src, P_src, flash, scp],\n connections=[\n Connection(F_src, flash), # F -> port 0\n Connection(z_src, flash[1]), # z_1 -> port 1\n Connection(T_src, flash[2]), # T -> port 2\n Connection(P_src, flash[3]), # P -> port 3\n Connection(flash, scp), # V_rate -> scope 0\n Connection(flash[1], scp[1]), # L_rate -> scope 1\n Connection(flash[2], scp[2]), # y_1 -> scope 2\n Connection(flash[3], scp[3]), # x_1 -> scope 3\n ],\n dt=0.5,\n)\n\nsim.run(120)" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "time, signals = scp.read()\nT_sweep = 340.0 + time * 0.5\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n\nax1.plot(T_sweep, signals[0], label=\"Vapor rate\")\nax1.plot(T_sweep, signals[1], label=\"Liquid rate\")\nax1.set_xlabel(\"Temperature [K]\")\nax1.set_ylabel(\"Flow rate [mol/s]\")\nax1.set_title(\"Flash Drum Flow Rates\")\nax1.legend()\nax1.grid(True, alpha=0.3)\n\nax2.plot(T_sweep, signals[2], label=r\"$y_1$ (vapor)\")\nax2.plot(T_sweep, signals[3], label=r\"$x_1$ (liquid)\")\nax2.axhline(0.5, color=\"gray\", ls=\"--\", alpha=0.4, label=\"Feed\")\nax2.set_xlabel(\"Temperature [K]\")\nax2.set_ylabel(\"Mole fraction (benzene)\")\nax2.set_title(\"VLE Compositions\")\nax2.legend()\nax2.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As temperature increases, more liquid vaporizes (higher vapor rate). The vapor is enriched in the lighter component (benzene), while the liquid becomes richer in toluene — the classic VLE separation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Part 2: Distillation Column (5 Trays)\n", + "\n", + "A distillation column is built by wiring multiple `DistillationTray` blocks in series. Each tray enforces vapor-liquid equilibrium with a constant relative volatility $\\alpha$:\n", + "\n", + "$$y = \\frac{\\alpha\\, x}{1 + (\\alpha - 1)\\, x}$$\n", + "\n", + "Under constant molar overflow (CMO), liquid flows down and vapor flows up with constant rates $L$ and $V$. We feed the column at the middle tray." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Column parameters\n", + "N_trays = 5\n", + "alpha = 2.5 # relative volatility\n", + "L = 5.0 # liquid flow rate [mol/s]\n", + "V = 5.0 # vapor flow rate [mol/s]\n", + "x_feed = 0.5 # feed composition (light component)\n", + "\n", + "# Create tray blocks (all start at x = 0.5)\n", + "trays = [DistillationTray(M=1.0, alpha=alpha, x0=0.5) for _ in range(N_trays)]\n", + "\n", + "# Liquid feed from condenser (enters tray 0 from above)\n", + "L_src = Source(func=lambda t: L)\n", + "x_top = Source(func=lambda t: 0.95) # reflux composition (nearly pure light)\n", + "\n", + "# Vapor feed from reboiler (enters tray N-1 from below)\n", + "V_src = Source(func=lambda t: V)\n", + "y_bot = Source(func=lambda t: 0.05) # reboiler vapor (nearly pure heavy)\n", + "\n", + "# Record composition on each tray\n", + "scp = Scope(labels=[f\"Tray {i+1}\" for i in range(N_trays)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Wire the trays: liquid cascades downward (tray $i$ liquid output $\\to$ tray $i+1$ liquid input), vapor rises upward (tray $i$ vapor output $\\to$ tray $i-1$ vapor input). External feeds enter the top and bottom." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "connections = []\n", + "\n", + "# Top tray (0): liquid from reflux\n", + "connections.append(Connection(L_src, trays[0])) # L_in -> port 0\n", + "connections.append(Connection(x_top, trays[0][1])) # x_in -> port 1\n", + "\n", + "# Bottom tray (N-1): vapor from reboiler\n", + "connections.append(Connection(V_src, trays[-1][2])) # V_in -> port 2\n", + "connections.append(Connection(y_bot, trays[-1][3])) # y_in -> port 3\n", + "\n", + "# Inter-tray connections\n", + "for i in range(N_trays - 1):\n", + " # Liquid flows down: tray i -> tray i+1\n", + " connections.append(Connection(trays[i], trays[i+1])) # L_out -> L_in\n", + " connections.append(Connection(trays[i][1], trays[i+1][1])) # x_out -> x_in\n", + "\n", + " # Vapor flows up: tray i+1 -> tray i\n", + " connections.append(Connection(trays[i+1][2], trays[i][2])) # V_out -> V_in\n", + " connections.append(Connection(trays[i+1][3], trays[i][3])) # y_out -> y_in\n", + "\n", + "# Connect each tray's liquid composition to scope\n", + "for i, tray in enumerate(trays):\n", + " connections.append(Connection(tray[1], scp[i])) # x_out -> scope\n", + "\n", + "sim = Simulation(\n", + " blocks=[L_src, x_top, V_src, y_bot, *trays, scp],\n", + " connections=connections,\n", + " dt=0.05,\n", + ")\n", + "\n", + "sim.run(30)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time, signals = scp.read()\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + "\n", + "# Dynamic tray composition evolution\n", + "for i in range(N_trays):\n", + " ax1.plot(time, signals[i], label=f\"Tray {i+1}\")\n", + "ax1.set_xlabel(\"Time [s]\")\n", + "ax1.set_ylabel(r\"$x$ (light component)\")\n", + "ax1.set_title(\"Tray Compositions Over Time\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Steady-state composition profile\n", + "x_profile = [tray.engine.state[0] for tray in trays]\n", + "tray_nums = list(range(1, N_trays + 1))\n", + "\n", + "ax2.plot(tray_nums, x_profile, \"o-\", markersize=8)\n", + "ax2.set_xlabel(\"Tray number (top to bottom)\")\n", + "ax2.set_ylabel(r\"$x$ (light component)\")\n", + "ax2.set_title(\"Composition Profile (Steady State)\")\n", + "ax2.grid(True, alpha=0.3)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The column separates the light and heavy components across its trays. The top tray is enriched in the light component (high $x$) while the bottom tray is depleted. The steady-state composition profile shows the characteristic staircase that a McCabe-Thiele diagram would predict for this relative volatility." + ] + } + ], + "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/heat_exchanger.ipynb b/docs/source/examples/heat_exchanger.ipynb new file mode 100644 index 0000000..7cd9287 --- /dev/null +++ b/docs/source/examples/heat_exchanger.ipynb @@ -0,0 +1,126 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Counter-Current Heat Exchanger\n", + "\n", + "Simulating the startup transient of a counter-current shell-and-tube heat exchanger, showing how temperature profiles develop along the exchanger length." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The exchanger is discretized into $N$ cells along its length. In each cell, the hot and cold streams exchange heat proportional to the local temperature difference:\n", + "\n", + "$$\\frac{dT_{h,i}}{dt} = \\frac{F_h}{V_\\text{cell,h}} (T_{h,i-1} - T_{h,i}) - \\frac{UA_\\text{cell}}{\\rho_h C_{p,h} V_\\text{cell,h}} (T_{h,i} - T_{c,i})$$\n", + "\n", + "$$\\frac{dT_{c,i}}{dt} = \\frac{F_c}{V_\\text{cell,c}} (T_{c,i+1} - T_{c,i}) + \\frac{UA_\\text{cell}}{\\rho_c C_{p,c} V_\\text{cell,c}} (T_{h,i} - T_{c,i})$$\n", + "\n", + "The hot stream flows from cell 1 to $N$, the cold stream flows counter-currently from cell $N$ to 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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.process import HeatExchanger" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "Configure a heat exchanger with 10 cells. Both sides start at 300 K (cold), then hot fluid at 370 K enters. This lets us watch the temperature profile build up from the initial uniform state. The high $UA$ relative to the flow rate ensures significant heat transfer." + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "hx = HeatExchanger(\n N_cells=10, # spatial discretization\n F_h=0.002, # hot flow rate [m³/s]\n F_c=0.002, # cold flow rate [m³/s]\n V_h=0.1, # hot-side volume [m³]\n V_c=0.1, # cold-side volume [m³]\n UA=5000.0, # overall heat transfer [W/K]\n rho_h=1000.0, # hot-side density [kg/m³]\n Cp_h=4184.0, # hot-side heat capacity [J/(kg·K)]\n rho_c=1000.0, # cold-side density [kg/m³]\n Cp_c=4184.0, # cold-side heat capacity [J/(kg·K)]\n T_h0=300.0, # initial hot-side temp [K]\n T_c0=300.0, # initial cold-side temp [K]\n)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Feed hot water at 370 K and cold water at 290 K. Record the outlet temperatures over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "T_h_in = Source(func=lambda t: 370.0) # hot inlet [K]\nT_c_in = Source(func=lambda t: 290.0) # cold inlet [K]\n\nscp = Scope(labels=[\"T_h_out [K]\", \"T_c_out [K]\"])\n\nsim = Simulation(\n blocks=[T_h_in, T_c_in, hx, scp],\n connections=[\n Connection(T_h_in, hx), # T_h_in -> port 0\n Connection(T_c_in, hx[1]), # T_c_in -> port 1\n Connection(hx, scp), # T_h_out -> scope port 0\n Connection(hx[1], scp[1]), # T_c_out -> scope port 1\n ],\n dt=0.02,\n)\n\nsim.run(200)" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time, signals = scp.read()\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + "\n", + "# Outlet temperatures over time\n", + "ax1.plot(time, signals[0], label=\"Hot outlet\")\n", + "ax1.plot(time, signals[1], label=\"Cold outlet\")\n", + "ax1.axhline(370, color=\"gray\", ls=\"--\", alpha=0.4, label=\"Hot inlet\")\n", + "ax1.axhline(290, color=\"gray\", ls=\"-.\", alpha=0.4, label=\"Cold inlet\")\n", + "ax1.set_xlabel(\"Time [s]\")\n", + "ax1.set_ylabel(\"Temperature [K]\")\n", + "ax1.set_title(\"Outlet Temperatures\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Steady-state temperature profile along the exchanger\n", + "state = hx.engine.x0\n", + "T_h_profile = state[0::2]\n", + "T_c_profile = state[1::2]\n", + "cells = np.arange(1, len(T_h_profile) + 1)\n", + "\n", + "ax2.plot(cells, T_h_profile, \"o-\", label=\"Hot stream\")\n", + "ax2.plot(cells, T_c_profile, \"s-\", label=\"Cold stream\")\n", + "ax2.set_xlabel(\"Cell number\")\n", + "ax2.set_ylabel(\"Temperature [K]\")\n", + "ax2.set_title(\"Temperature Profile (Steady State)\")\n", + "ax2.legend()\n", + "ax2.grid(True, alpha=0.3)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At startup, both sides are cold. The hot fluid quickly warms the exchanger from the inlet side while the cold fluid absorbs heat. At steady state, the counter-current arrangement creates the characteristic temperature cross — the cold outlet approaches the hot inlet and vice versa, achieving high thermal efficiency." + ] + } + ], + "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/src/pathsim_chem/process/distillation.py b/src/pathsim_chem/process/distillation.py index 1656f06..5371ee7 100644 --- a/src/pathsim_chem/process/distillation.py +++ b/src/pathsim_chem/process/distillation.py @@ -82,8 +82,18 @@ def __init__(self, M=1.0, alpha=2.5, x0=0.5): def _vle(x_val): return self.alpha * x_val / (1.0 + (self.alpha - 1.0) * x_val) + # ensure u has expected 4 elements (handles framework probing) + def _pad_u(u): + u = np.atleast_1d(u) + if len(u) < 4: + padded = np.zeros(4) + padded[:len(u)] = u + return padded + return u + # rhs of tray ode def _fn_d(x, u, t): + u = _pad_u(u) L_in, x_in, V_in, y_in = u x_tray = x[0] @@ -95,6 +105,7 @@ def _fn_d(x, u, t): # jacobian wrt x def _jc_d(x, u, t): + u = _pad_u(u) L_in, x_in, V_in, y_in = u x_tray = x[0] @@ -106,6 +117,7 @@ def _jc_d(x, u, t): # output: L_out, x, V_out, y (CMO: pass-through flows) def _fn_a(x, u, t): + u = _pad_u(u) L_in, x_in, V_in, y_in = u x_tray = x[0] y_tray = _vle(x_tray) diff --git a/src/pathsim_chem/process/flash_drum.py b/src/pathsim_chem/process/flash_drum.py index 33c5f43..29f71a0 100644 --- a/src/pathsim_chem/process/flash_drum.py +++ b/src/pathsim_chem/process/flash_drum.py @@ -87,8 +87,9 @@ def __init__(self, holdup=100.0, self.holdup = holdup - # default Antoine parameters: benzene / toluene (Pa, K) - self.antoine_A = np.array(antoine_A) if antoine_A is not None else np.array([15.9008, 16.0137]) + # default Antoine parameters: benzene / toluene (ln(Pa), K) + # Converted from log10/mmHg: A_ln = ln(10)*A_log10 + ln(133.322) + self.antoine_A = np.array(antoine_A) if antoine_A is not None else np.array([20.7936, 20.9064]) self.antoine_B = np.array(antoine_B) if antoine_B is not None else np.array([2788.51, 3096.52]) self.antoine_C = np.array(antoine_C) if antoine_C is not None else np.array([-52.36, -53.67]) @@ -107,16 +108,23 @@ def _solve_vle(z, T, P): P_sat = np.exp(self.antoine_A - self.antoine_B / (T + self.antoine_C)) K = P_sat / P - d1 = K[0] - 1.0 - d2 = K[1] - 1.0 + # check bubble/dew point conditions + bubble = np.sum(z * K) # if < 1, subcooled (all liquid) + dew = np.sum(z / K) if np.all(K > 1e-12) else np.inf # if < 1, superheated (all vapor) - if abs(d1) < 1e-12 and abs(d2) < 1e-12: + if bubble <= 1.0: + # subcooled: all liquid beta = 0.0 + elif dew <= 1.0: + # superheated: all vapor + beta = 1.0 else: + # two-phase: solve Rachford-Rice + d1 = K[0] - 1.0 + d2 = K[1] - 1.0 den = d1 * d2 beta = -(z[0] * d1 + z[1] * d2) / den if abs(den) > 1e-12 else 0.0 - - beta = np.clip(beta, 0.0, 1.0) + beta = np.clip(beta, 0.0, 1.0) y = z * K / (1.0 + beta * (K - 1.0)) x_eq = z / (1.0 + beta * (K - 1.0)) @@ -130,8 +138,18 @@ def _solve_vle(z, T, P): return beta, y, x_eq + # ensure u has expected 4 elements (handles framework probing) + def _pad_u(u): + u = np.atleast_1d(u) + if len(u) < 4: + padded = np.zeros(4) + padded[:len(u)] = u + return padded + return u + # rhs of flash drum ode def _fn_d(x, u, t): + u = _pad_u(u) F_in, z_1, T, P = u z = np.array([z_1, 1.0 - z_1]) @@ -144,6 +162,7 @@ def _fn_d(x, u, t): # algebraic output def _fn_a(x, u, t): + u = _pad_u(u) F_in, z_1, T, P = u z = np.array([z_1, 1.0 - z_1]) diff --git a/src/pathsim_chem/process/heat_exchanger.py b/src/pathsim_chem/process/heat_exchanger.py index 2e5465c..4c15cfb 100644 --- a/src/pathsim_chem/process/heat_exchanger.py +++ b/src/pathsim_chem/process/heat_exchanger.py @@ -112,8 +112,18 @@ def __init__(self, N_cells=5, F_h=0.1, F_c=0.1, V_h=0.5, V_c=0.5, x0[0::2] = T_h0 x0[1::2] = T_c0 + # ensure u has expected 2 elements (handles framework probing) + def _pad_u(u): + u = np.atleast_1d(u) + if len(u) < 2: + padded = np.zeros(2) + padded[:len(u)] = u + return padded + return u + # rhs of heat exchanger ode (vectorized) def _fn_d(x, u, t): + u = _pad_u(u) T_h_in, T_c_in = u N = self.N_cells