Skip to content

Commit 199625a

Browse files
milanofthemilanofthe
committed
Add PFR block with discretized Arrhenius kinetics and energy balance
1 parent 8cb3016 commit 199625a

3 files changed

Lines changed: 370 additions & 0 deletions

File tree

‎src/pathsim_chem/process/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@
99
from .flash_drum import *
1010
from .distillation import *
1111
from .multicomponent_flash import *
12+
from .pfr import *

‎src/pathsim_chem/process/pfr.py‎

Lines changed: 209 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,209 @@
1+
#########################################################################################
2+
##
3+
## Plug Flow Reactor (PFR) Block
4+
##
5+
#########################################################################################
6+
7+
# IMPORTS ===============================================================================
8+
9+
import numpy as np
10+
11+
from pathsim.blocks.dynsys import DynamicalSystem
12+
13+
# CONSTANTS =============================================================================
14+
15+
R_GAS = 8.314 # universal gas constant [J/(mol·K)]
16+
17+
# BLOCKS ================================================================================
18+
19+
class PFR(DynamicalSystem):
20+
"""Plug flow reactor with Arrhenius kinetics and energy balance.
21+
22+
Discretized tubular reactor divided into N cells along its length.
23+
Each cell has concentration and temperature states with nth-order
24+
kinetics and an energy balance including heat of reaction.
25+
26+
Mathematical Formulation
27+
-------------------------
28+
For each cell :math:`i = 1, \\ldots, N`:
29+
30+
.. math::
31+
32+
\\frac{dC_i}{dt} = \\frac{F}{V_{cell}} (C_{i-1} - C_i) - k(T_i) \\, C_i^n
33+
34+
.. math::
35+
36+
\\frac{dT_i}{dt} = \\frac{F}{V_{cell}} (T_{i-1} - T_i)
37+
+ \\frac{(-\\Delta H_{rxn})}{\\rho \\, C_p} \\, k(T_i) \\, C_i^n
38+
39+
where the Arrhenius rate constant is:
40+
41+
.. math::
42+
43+
k(T) = k_0 \\, \\exp\\!\\left(-\\frac{E_a}{R \\, T}\\right)
44+
45+
The state vector is ordered as
46+
:math:`[C_1, T_1, C_2, T_2, \\ldots, C_N, T_N]`.
47+
48+
Parameters
49+
----------
50+
N_cells : int
51+
Number of discretization cells [-].
52+
V : float
53+
Total reactor volume [m³].
54+
F : float
55+
Volumetric flow rate [m³/s].
56+
k0 : float
57+
Pre-exponential Arrhenius factor [1/s for n=1].
58+
Ea : float
59+
Activation energy [J/mol].
60+
n : float
61+
Reaction order [-].
62+
dH_rxn : float
63+
Heat of reaction [J/mol]. Negative for exothermic.
64+
rho : float
65+
Fluid density [kg/m³].
66+
Cp : float
67+
Fluid heat capacity [J/(kg·K)].
68+
C0 : float
69+
Initial concentration [mol/m³].
70+
T0 : float
71+
Initial temperature [K].
72+
"""
73+
74+
input_port_labels = {
75+
"C_in": 0,
76+
"T_in": 1,
77+
}
78+
79+
output_port_labels = {
80+
"C_out": 0,
81+
"T_out": 1,
82+
}
83+
84+
def __init__(self, N_cells=5, V=1.0, F=0.1, k0=1e6, Ea=50000.0, n=1.0,
85+
dH_rxn=-50000.0, rho=1000.0, Cp=4184.0,
86+
C0=0.0, T0=300.0):
87+
88+
# input validation
89+
if N_cells < 1:
90+
raise ValueError(f"'N_cells' must be >= 1 but is {N_cells}")
91+
if V <= 0:
92+
raise ValueError(f"'V' must be positive but is {V}")
93+
if F <= 0:
94+
raise ValueError(f"'F' must be positive but is {F}")
95+
if rho <= 0:
96+
raise ValueError(f"'rho' must be positive but is {rho}")
97+
if Cp <= 0:
98+
raise ValueError(f"'Cp' must be positive but is {Cp}")
99+
100+
# store parameters
101+
self.N_cells = int(N_cells)
102+
self.V = V
103+
self.F = F
104+
self.k0 = k0
105+
self.Ea = Ea
106+
self.n = n
107+
self.dH_rxn = dH_rxn
108+
self.rho = rho
109+
self.Cp = Cp
110+
111+
N = self.N_cells
112+
113+
# initial state: interleaved [C_1, T_1, C_2, T_2, ...]
114+
x0 = np.empty(2 * N)
115+
x0[0::2] = C0
116+
x0[1::2] = T0
117+
118+
# ensure u has expected 2 elements (handles framework probing)
119+
def _pad_u(u):
120+
u = np.atleast_1d(u)
121+
if len(u) < 2:
122+
padded = np.zeros(2)
123+
padded[:len(u)] = u
124+
return padded
125+
return u
126+
127+
# rhs of PFR ode (vectorized)
128+
def _fn_d(x, u, t):
129+
u = _pad_u(u)
130+
C_in, T_in = u
131+
N = self.N_cells
132+
133+
V_cell = self.V / N
134+
f_flow = self.F / V_cell
135+
rcp = (-self.dH_rxn) / (self.rho * self.Cp)
136+
137+
C = x[0::2]
138+
T = x[1::2]
139+
140+
# upstream values with boundary conditions
141+
C_prev = np.empty(N)
142+
C_prev[0] = C_in
143+
C_prev[1:] = C[:-1]
144+
145+
T_prev = np.empty(N)
146+
T_prev[0] = T_in
147+
T_prev[1:] = T[:-1]
148+
149+
# Arrhenius rate per cell
150+
k = self.k0 * np.exp(-self.Ea / (R_GAS * T))
151+
r = k * np.abs(C)**self.n # abs for numerical safety
152+
153+
dx = np.empty(2 * N)
154+
dx[0::2] = f_flow * (C_prev - C) - r
155+
dx[1::2] = f_flow * (T_prev - T) + rcp * r
156+
157+
return dx
158+
159+
# analytical jacobian (block-tridiagonal structure)
160+
def _jc_d(x, u, t):
161+
N = self.N_cells
162+
163+
V_cell = self.V / N
164+
f_flow = self.F / V_cell
165+
rcp = (-self.dH_rxn) / (self.rho * self.Cp)
166+
167+
C = x[0::2]
168+
T = x[1::2]
169+
170+
k = self.k0 * np.exp(-self.Ea / (R_GAS * T))
171+
dk_dT = k * self.Ea / (R_GAS * T**2)
172+
173+
dim = 2 * N
174+
J = np.zeros((dim, dim))
175+
176+
for i in range(N):
177+
ci = 2 * i # concentration index
178+
ti = 2 * i + 1 # temperature index
179+
180+
C_i = max(abs(C[i]), 1e-30)
181+
dr_dC = k[i] * self.n * C_i**(self.n - 1) if C_i > 0 else 0.0
182+
dr_dT = dk_dT[i] * C_i**self.n
183+
184+
# dC_i/dC_i, dC_i/dT_i
185+
J[ci, ci] = -f_flow - dr_dC
186+
J[ci, ti] = -dr_dT
187+
188+
# dT_i/dC_i, dT_i/dT_i
189+
J[ti, ci] = rcp * dr_dC
190+
J[ti, ti] = -f_flow + rcp * dr_dT
191+
192+
# upstream coupling: dC_i/dC_{i-1}, dT_i/dT_{i-1}
193+
if i > 0:
194+
J[ci, 2*(i-1)] = f_flow
195+
J[ti, 2*(i-1) + 1] = f_flow
196+
197+
return J
198+
199+
# output: last cell values
200+
def _fn_a(x, u, t):
201+
N = self.N_cells
202+
return np.array([x[2*(N-1)], x[2*(N-1) + 1]])
203+
204+
super().__init__(
205+
func_dyn=_fn_d,
206+
jac_dyn=_jc_d,
207+
func_alg=_fn_a,
208+
initial_value=x0,
209+
)

‎tests/process/test_pfr.py‎

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
########################################################################################
2+
##
3+
## TESTS FOR
4+
## 'process.pfr.py'
5+
##
6+
########################################################################################
7+
8+
# IMPORTS ==============================================================================
9+
10+
import unittest
11+
import numpy as np
12+
13+
from pathsim_chem.process import PFR
14+
15+
from pathsim.solvers import EUF
16+
17+
18+
# TESTS ================================================================================
19+
20+
class TestPFR(unittest.TestCase):
21+
"""Test the plug flow reactor block."""
22+
23+
def test_init_default(self):
24+
"""Test default initialization."""
25+
P = PFR()
26+
self.assertEqual(P.N_cells, 5)
27+
self.assertEqual(P.V, 1.0)
28+
self.assertEqual(P.F, 0.1)
29+
self.assertEqual(P.k0, 1e6)
30+
self.assertEqual(P.Ea, 50000.0)
31+
self.assertEqual(P.n, 1.0)
32+
33+
def test_init_custom(self):
34+
"""Test custom initialization."""
35+
P = PFR(N_cells=10, V=2.0, F=0.5, k0=1e4, Ea=40000.0, n=2.0,
36+
dH_rxn=-30000.0, rho=800.0, Cp=3000.0, C0=1.5, T0=350.0)
37+
self.assertEqual(P.N_cells, 10)
38+
self.assertEqual(P.V, 2.0)
39+
40+
P.set_solver(EUF, parent=None)
41+
state = P.engine.initial_value
42+
self.assertEqual(len(state), 20) # 2 * N_cells
43+
self.assertTrue(np.all(state[0::2] == 1.5)) # concentrations
44+
self.assertTrue(np.all(state[1::2] == 350.0)) # temperatures
45+
46+
def test_init_validation(self):
47+
"""Test input validation."""
48+
with self.assertRaises(ValueError):
49+
PFR(N_cells=0)
50+
with self.assertRaises(ValueError):
51+
PFR(V=-1)
52+
with self.assertRaises(ValueError):
53+
PFR(F=0)
54+
with self.assertRaises(ValueError):
55+
PFR(rho=-1)
56+
with self.assertRaises(ValueError):
57+
PFR(Cp=0)
58+
59+
def test_port_labels(self):
60+
"""Test port label definitions."""
61+
self.assertEqual(PFR.input_port_labels["C_in"], 0)
62+
self.assertEqual(PFR.input_port_labels["T_in"], 1)
63+
self.assertEqual(PFR.output_port_labels["C_out"], 0)
64+
self.assertEqual(PFR.output_port_labels["T_out"], 1)
65+
66+
def test_state_size(self):
67+
"""Test that state vector has correct size."""
68+
for N in [1, 3, 10]:
69+
P = PFR(N_cells=N)
70+
P.set_solver(EUF, parent=None)
71+
self.assertEqual(len(P.engine.initial_value), 2 * N)
72+
73+
def test_output_initial(self):
74+
"""Test outputs at initial state."""
75+
P = PFR(N_cells=3, C0=2.0, T0=350.0)
76+
P.set_solver(EUF, parent=None)
77+
P.update(None)
78+
79+
# Last cell values
80+
self.assertAlmostEqual(P.outputs[0], 2.0) # C_out
81+
self.assertAlmostEqual(P.outputs[1], 350.0) # T_out
82+
83+
def test_no_reaction(self):
84+
"""With k0=0, no reaction occurs. At steady state C_out = C_in."""
85+
P = PFR(N_cells=3, V=1.0, F=1.0, k0=0.0, Ea=0.0, n=1.0,
86+
dH_rxn=0.0, C0=1.0, T0=350.0)
87+
P.set_solver(EUF, parent=None)
88+
89+
P.inputs[0] = 1.0 # C_in
90+
P.inputs[1] = 350.0 # T_in
91+
92+
# When C=C_in and T=T_in everywhere, with no reaction, derivatives = 0
93+
x = P.engine.get()
94+
u = np.array([1.0, 350.0])
95+
dx = P.op_dyn(x, u, 0)
96+
97+
self.assertTrue(np.allclose(dx, 0.0, atol=1e-10))
98+
99+
def test_reaction_direction(self):
100+
"""With reaction, concentration should decrease along the reactor."""
101+
R_GAS = 8.314
102+
T = 350.0
103+
k0, Ea = 1e6, 50000.0
104+
k = k0 * np.exp(-Ea / (R_GAS * T))
105+
106+
P = PFR(N_cells=1, V=1.0, F=1.0, k0=k0, Ea=Ea, n=1.0,
107+
dH_rxn=0.0, C0=1.0, T0=T)
108+
P.set_solver(EUF, parent=None)
109+
110+
# Set inlet = current state, so flow terms cancel
111+
P.inputs[0] = 1.0 # C_in
112+
P.inputs[1] = T # T_in
113+
114+
x = np.array([1.0, T])
115+
u = np.array([1.0, T])
116+
dx = P.op_dyn(x, u, 0)
117+
118+
# dC/dt = 0 (flow) - k*C = -k
119+
self.assertAlmostEqual(dx[0], -k, places=5)
120+
# dT/dt = 0 since dH_rxn=0
121+
self.assertAlmostEqual(dx[1], 0.0, places=5)
122+
123+
def test_energy_conservation_exothermic(self):
124+
"""For exothermic reaction (dH_rxn < 0), temperature should increase."""
125+
P = PFR(N_cells=1, V=1.0, F=1.0, k0=1e3, Ea=20000.0, n=1.0,
126+
dH_rxn=-50000.0, rho=1000.0, Cp=4184.0,
127+
C0=2.0, T0=350.0)
128+
P.set_solver(EUF, parent=None)
129+
130+
P.inputs[0] = 2.0
131+
P.inputs[1] = 350.0
132+
133+
x = np.array([2.0, 350.0])
134+
u = np.array([2.0, 350.0])
135+
dx = P.op_dyn(x, u, 0)
136+
137+
# Concentration should decrease (reaction consuming)
138+
self.assertLess(dx[0], 0)
139+
# Temperature should increase (exothermic)
140+
self.assertGreater(dx[1], 0)
141+
142+
def test_jacobian_stability(self):
143+
"""Diagonal of Jacobian should be negative (stable)."""
144+
P = PFR(N_cells=2, V=1.0, F=0.1, k0=1e3, Ea=30000.0, n=1.0,
145+
dH_rxn=0.0, C0=1.0, T0=350.0)
146+
P.set_solver(EUF, parent=None)
147+
148+
x = P.engine.get()
149+
u = np.array([1.0, 350.0])
150+
J = P.op_dyn.jac_x(x, u, 0)
151+
152+
# All diagonal elements should be negative
153+
for i in range(len(x)):
154+
self.assertLess(J[i, i], 0)
155+
156+
157+
# RUN TESTS LOCALLY ====================================================================
158+
159+
if __name__ == '__main__':
160+
unittest.main(verbosity=2)

0 commit comments

Comments
 (0)