Dynamic Operability Analysis of a Membrane Reactor for Methane Aromatization and Hydrogen Production#
Author: Victor Alves - Carnegie Mellon University
This case study consists of a membrane reactor for direct methane aromatization (DMA-MR) that allows hydrogen and benzene production from natural gas. This single-unit operation is capable of performing reaction and separation simultaneously, achieving higher conversions due to Le Chatelier’s principle. The steady-state operability of the DMA-MR has been studied extensively. In practice, though, a modular unit runs continuously and its feed composition drifts over time. A dynamic operability analysis [10] answers a question the steady-state picture cannot: starting from the current operating point, can we move the inputs fast enough to steer the outputs where we want them, and to absorb feed disturbances, within a useful time window?
In this case study, some of opyrability’s features are highlighted, such as:
Building the dynamic achievable output funnel (\(AOS^k\)): the set of outputs the reactor can reach within \(k\) time steps, for a nonlinear model with many internal states (160 here)
Computing the dynamic Operability Index (dOI) at each time step (how much of the desired output set is reachable so far) and showing it directly on the funnel
Checking robustness: which outputs stay reachable across a range of methane feed disturbances
The schematic below [3] depicts the membrane reactor in a nutshell:
The two-step reaction mechanism for the nonoxidative conversion of methane is:
For details on the first-principles model steady-state equations, see the Membrane Reactor example in this gallery. Here, the dynamic model from [10] is used: the reactor length is discretized into \(N = 20\) elements and the state is the spatial profile of molar holdups (4 tube and 4 shell species per element, 160 states in total), integrated in time as a system of ODEs. The model is defined in full below so this notebook is self-contained (it reproduces the same model that ships in opyrability’s test suite, tests/dma_mr.py), providing three functions: one that computes the reactor’s starting steady state, one that advances the reactor by a single time step, and one that reads the outputs off the reactor state.
Variables:
Input (AIS) |
Output (AOS) |
Disturbance (EDS) |
|---|---|---|
Tube flowrate [dm³/h] |
Benzene mole fraction in tube outlet [%] |
Methane feed mole fraction [%] |
Shell flowrate [dm³/h] |
Hydrogen mole fraction in shell outlet [%] |
In this case study, the dynamic achievable output funnel will be mapped from the nominal steady state, and the dOI will be evaluated against a desired output set at each time step. Then, the funnel will be re-evaluated under two methane feed disturbance scenarios to quantify the outputs that are achievable regardless of the disturbance realization.
The dynamic model equations#
The reactor is described by transient mole balances on the molar holdups of each species, in both the tube and the shell, resolved along the reactor length. The two reactions above proceed at equilibrium-limited rates
where \(v_b\) is the bed voidage, \(\eta\) the catalyst effectiveness, \(k_1,k_2\) the rate constants and \(K_1,K_2\) the equilibrium constants. Removing hydrogen through the membrane shifts both equilibria forward (Le Chatelier).
Membrane permeation: Each species \(j\) permeates from tube to shell with a quarter-power (Sieverts-type) driving force,
where \(Q\) is the membrane permeance and \(\alpha_j\) the selectivity (\(\alpha = 1\) for \(\mathrm{H_2}\), which permeates freely, and \(\alpha = 100\) for the other species).
Spatial discretization: The reactor of length \(L = N\,\Delta z\) is divided into \(N = 20\) well-mixed finite-volume elements (a tanks-in-series approximation of the plug-flow profile). Element \(i\) holds the tube amounts \(M_{t,j,i}\) and shell amounts \(M_{s,j,i}\) (four species each, so \(8N = 160\) states). Neighbouring elements are coupled by convection: element \(i\) receives the molar outflow of element \(i-1\). Because the reactor is co-current, the tube feed and the shell sweep both enter at element \(1\) (\(z=0\)) and flow toward element \(N\) (\(z=L\)).
Dynamic mole balances: For each element \(i\) and species \(j\),
with the tube inlet \(F_{t,j,0} = Q_t\,C_t\,y_{j}^{\text{feed}}\) and the net reaction-generation coefficients
The convective outflows treat the non-reacting, non-permeating inert (tube) and sweep gas (shell) as a tie component (a stream whose flow is fixed, used to back out the others), for example \(F_{t,j,i} = -F_{t,\mathrm{inert}}\,M_{t,j,i} \big/ \left(\sum_j M_{t,j,i} - C_t V_t\right)\). All balances are scaled by a time factor so that integrating over one unit of time advances the reactor by one minute.
How the system is solved:
Initial steady state: The steady-state spatial profile \(dF/dz\) is integrated once along \(z\) for a good guess, the molar flows are converted to holdups, and the transient model is then relaxed to steady state (the reactor settles to a steady state on its own, so simply marching forward in time converges to it).
Funnel: Starting from that steady state, the \(160\)-state stiff ODE system is integrated one minute per step (with
scipy.integrate.solve_ivp, LSODA method) while the inlet flowrates are held at each candidate value. Stacking the reachable outputs across the time steps builds the dynamic funnel.
Outputs: The two monitored outputs are taken at the reactor end (\(z=L\), element \(N\)):
Let’s start by importing the packages and defining the DMA-MR dynamic model:
from __future__ import annotations
import numpy as np
import scipy.integrate as spint
from numpy import pi as pi
import matplotlib.pyplot as plt
# Plotly renderer that stays interactive in the compiled
# documentation website (loads plotly.js via require.js).
import plotly.io as pio
pio.renderers.default = 'notebook_connected'
# --------------------------------------------------------------------------- #
# Co-current dynamic DMA-MR model
#
# Co-current shell-and-tube membrane reactor: the tube feed and the shell sweep
# gas both enter at z = 0 and flow the same way to z = L. The state is the
# spatial profile of molar holdups, n_elements elements times 8 species (4 tube
# and 4 shell), stored as one flat vector. Units: dm, dm^3/h, mol.
# --------------------------------------------------------------------------- #
# Model parameters (defaults reproduce the published case study).
P = 1.0 # Operating pressure [atm].
Ct = 0.010387973185676 # Total concentration P/(R T) [mol/dm^3].
Dt = 1.5 # Tube diameter [dm].
Vt = 4.417864669110937 # Per-element tube volume [dm^3].
Vs = 1.943860454408813 # Per-element shell volume [dm^3].
dz = 2.5 # Element length [dm]; reactor length L = N * dz.
vb = 0.5 # Catalyst bed voidage.
eff = 0.9 # Catalyst effectiveness factor.
k1 = 144.0 # Reaction 1 forward rate constant [1/h].
k2 = 15120.0 # Reaction 2 forward rate constant [1/h].
k1eq = 6.25e-06 # Reaction 1 equilibrium constant.
k2eq = 74.494 # Reaction 2 equilibrium constant.
permeance = 0.01 # Membrane permeance (applied as permeance * 36).
alpha = 100.0 # Membrane selectivity (H2 over the other species).
N = 20 # Number of spatial elements.
timescale = 1.0 / 60.0 # Integrating the model over [0, 1] is one minute.
feed_CH4 = 0.942 # Nominal CH4 feed mole fraction.
feed_C2H4 = 0.023 # Ethylene feed mole fraction (kept > 0).
feed_inert = 0.035 # Inert (N2) feed mole fraction.
def feed_fractions(d: float | None) -> tuple[float, float, float]:
"""
Map a disturbance ``d`` (methane feed mole percent) to the feed mole
fractions (CH4, C2H4, inert).
Ethylene is held at its nominal value and the inert takes up the balance.
``d`` may be a fraction (<= 1) or a percent (> 1); ``None`` returns the
nominal feed.
"""
if d is None:
return feed_CH4, feed_C2H4, feed_inert
yCH4 = d / 100.0 if d > 1.0 else float(d)
return yCH4, feed_C2H4, 1.0 - yCH4 - feed_C2H4
def steady_state_rhs(z: float, F: np.ndarray, Qt: float, Qs: float,
yCH4: float, yC2H4: float, yinert: float) -> list:
"""
Spatial steady-state molar-flow derivative dF/dz, integrated along the
reactor length to build the initial holdup profile.
``F`` holds the eight molar flows ``[tube CH4, C2H4, H2, C6H6, shell CH4,
C2H4, H2, C6H6]`` [mol/h]; ``Qt``/``Qs`` are the tube/shell inlet flowrates
[dm^3/h] and ``yCH4``/``yC2H4``/``yinert`` the feed mole fractions.
"""
def quarter_root(pressure):
# Guarded fourth root: the membrane driving force is P ** 0.25, and the
# solver can momentarily probe slightly negative pressures, for which a
# plain power would return NaN. Clipping keeps the integration real.
return np.power(np.clip(pressure, 1e-12, None), 0.25)
# Unpack the eight molar flows (4 tube, 4 shell).
FtCH4, FtC2H4, FtH2, FtC6H6, FsCH4, FsC2H4, FsH2, FsC6H6 = F
# Tube cross-sectional area and total molar flows (the inert is the tube
# tie component, the sweep gas the shell tie component).
At = Dt ** 2 / 4 * pi
Ftinert = Qt * Ct * yinert
Fshell0 = Qs * Ct
Ftube = FtCH4 + FtC2H4 + FtH2 + FtC6H6 + Ftinert
Fshell = FsCH4 + FsC2H4 + FsH2 + FsC6H6 + Fshell0
# Tube concentrations and tube/shell partial pressures.
CtCH4 = FtCH4 * Ct / Ftube
CtC2H4 = FtC2H4 * Ct / Ftube
CtH2 = FtH2 * Ct / Ftube
CtC6H6 = FtC6H6 * Ct / Ftube
Pt = (P * FtCH4 / Ftube, P * FtC2H4 / Ftube,
P * FtH2 / Ftube, P * FtC6H6 / Ftube)
Ps = (P * FsCH4 / Fshell, P * FsC2H4 / Fshell,
P * FsH2 / Fshell, P * FsC6H6 / Fshell)
# Membrane permeation from tube to shell. H2 permeates freely; the other
# species are slowed by the membrane selectivity.
perm_other = permeance * 36 / alpha
perm_H2 = permeance * 36
T2S_CH4 = perm_other * pi * Dt * (quarter_root(Pt[0]) - quarter_root(Ps[0]))
T2S_C2H4 = perm_other * pi * Dt * (quarter_root(Pt[1]) - quarter_root(Ps[1]))
T2S_H2 = perm_H2 * pi * Dt * (quarter_root(Pt[2]) - quarter_root(Ps[2]))
T2S_C6H6 = perm_other * pi * Dt * (quarter_root(Pt[3]) - quarter_root(Ps[3]))
# Equilibrium-limited rates of the two reactions.
r1 = -(1 - vb) * eff * k1 * CtCH4 * (
(k1eq * CtCH4 ** 2) - CtC2H4 * CtH2 ** 2) / (k1eq * CtCH4 ** 2)
r2 = -(1 - vb) * eff * k2 * CtC2H4 * (
(k2eq * CtC2H4 ** 3) - CtC6H6 * CtH2 ** 3) / (k2eq * CtC2H4 ** 3)
# dF/dz: four tube species (reaction generation minus permeation out) then
# four shell species (permeation in).
return [r1 * At - T2S_CH4,
(-r1 / 2 + r2) * At - T2S_C2H4,
(-r1 / 2 - r2) * At - T2S_H2,
(-r2 / 3) * At - T2S_C6H6,
T2S_CH4, T2S_C2H4, T2S_H2, T2S_C6H6]
def state_space_rhs(t: float, M: np.ndarray, Qt: float, Qs: float,
yCH4: float, yC2H4: float, yinert: float) -> np.ndarray:
"""
Transient state-space derivative dM/dt for the molar holdups, vectorized
over the spatial elements.
Integrating this right-hand side over one unit of time advances the reactor
by one minute. ``M`` is the flat holdup state of length ``8 * N``.
"""
def quarter_root(pressure):
# Guarded fourth root (see steady_state_rhs); here applied to the whole
# (N, 4) partial-pressure arrays at once.
return np.power(np.clip(pressure, 1e-12, None), 0.25)
# Inlet molar flows of the tube feed and the shell sweep gas.
Ft0 = Qt * Ct
Fs0 = Qs * Ct
Ftinert = Ft0 * yinert
feed = np.array([Ft0 * yCH4, Ft0 * yC2H4, 0.0, 0.0])
# Reshape the flat state into (element, species); split tube and shell.
X = np.asarray(M).reshape(N, 8)
Mt, Ms = X[:, 0:4], X[:, 4:8]
# Convective molar outflow of each element, with the non-permeating,
# non-reacting inert/sweep acting as the tie component (fixed flow).
denom_t = Mt.sum(axis=1) - Ct * Vt
denom_s = Ms.sum(axis=1) - Ct * Vs
Ft = -(Ftinert * Mt) / denom_t[:, None]
Fs = -(Fs0 * Ms) / denom_s[:, None]
# Tube concentrations and tube/shell partial pressures.
Ct_sp = Mt / Vt
Pt = Mt * P / (Ct * Vt)
Ps = Ms * P / (Ct * Vs)
# Membrane permeation from tube to shell (H2 unscaled, others divided by
# the selectivity).
sel = np.array([alpha, alpha, 1.0, alpha])
coef = (permeance * 36.0 / sel) * pi * Dt * dz
T2S = coef * (quarter_root(Pt) - quarter_root(Ps))
# Equilibrium-limited reaction rates per element.
CtCH4, CtC2H4, CtH2, CtC6H6 = Ct_sp[:, 0], Ct_sp[:, 1], Ct_sp[:, 2], Ct_sp[:, 3]
r1 = -(1 - vb) * eff * k1 * CtCH4 * (
(k1eq * CtCH4 ** 2) - CtC2H4 * CtH2 ** 2) / (k1eq * CtCH4 ** 2)
r2 = -(1 - vb) * eff * k2 * CtC2H4 * (
(k2eq * CtC2H4 ** 3) - CtC6H6 * CtH2 ** 3) / (k2eq * CtC2H4 ** 3)
# Net reaction generation of each tube species, scaled by the element
# volume (stoichiometry of 2 CH4 -> C2H4 + 2 H2 and 3 C2H4 -> C6H6 + 3 H2).
gen = np.empty((N, 4))
gen[:, 0] = r1
gen[:, 1] = -r1 / 2 + r2
gen[:, 2] = -r1 - r2
gen[:, 3] = -r2 / 3
gen *= Vt
# Inflow to each element from its upstream neighbor (co-current): tube
# element 0 is fed by the inlet, shell element 0 has no upstream neighbor.
Ft_in = np.vstack([feed, Ft[:-1, :]])
Fs_in = np.vstack([np.zeros(4), Fs[:-1, :]])
# Holdup balances: tube (convection + reaction - permeation out) and shell
# (convection + permeation in), scaled to a per-minute time base.
dM = np.empty((N, 8))
dM[:, 0:4] = (Ft_in - Ft + gen - T2S) * timescale
dM[:, 4:8] = (Fs_in - Fs + T2S) * timescale
return dM.reshape(-1)
def dma_mr_cocurrent_outputs(M: np.ndarray) -> np.ndarray:
"""
Reactor outputs from the holdup state ``M``: the benzene mole percent at
the tube outlet and the H2 mole percent at the shell outlet, both at the
reactor end (the last spatial element of the flat state).
"""
M = np.asarray(M)
# The last element block is M[-8:]: its tube C6H6 is M[-5] and shell H2 is
# M[-2]. Convert those holdups to mole percent.
benzene_pct = 100.0 * M[-5] / (Ct * Vt)
hydrogen_pct = 100.0 * M[-2] / (Ct * Vs)
return np.array([benzene_pct, hydrogen_pct])
def dma_mr_cocurrent_x0(u: np.ndarray, d: float | None = None,
n_relax: int = 50) -> np.ndarray:
"""
Steady-state holdup state at inputs ``u = [tube flowrate, shell flowrate]``
(dm^3/h) and disturbance ``d`` (methane feed mole percent; ``None`` for the
nominal feed).
A spatial-profile guess is built by integrating the steady-state balances
along the reactor, then relaxed to steady state with the transient model
(the reactor is Lyapunov-stable, so time-marching settles to steady state).
Returns the flat holdup state of length ``8 * N``.
"""
yCH4, yC2H4, yinert = feed_fractions(d)
Qt, Qs = float(u[0]), float(u[1])
# Inlet molar flows and the spatial grid along the reactor length.
Ft0 = Qt * Ct
Ftinert = Ft0 * yinert
Fs0 = Qs * Ct
F0 = [Ft0 * yCH4, Ft0 * yC2H4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Lgrid = np.linspace(0, dz * N, N + 1)
# Integrate the steady-state molar balances along z for a profile guess.
sol = spint.solve_ivp(steady_state_rhs, [0, Lgrid[-1]], F0, t_eval=Lgrid,
args=(Qt, Qs, yCH4, yC2H4, yinert),
rtol=1e-7, atol=1e-7, method='LSODA')
Fss = sol.y.T
# Convert the molar-flow profile to holdups (mole fractions times the
# element holdup Ct * V), dropping the inlet node.
M = np.zeros((Fss.shape[0], 8))
st = Fss[:, 0:4].sum(axis=1) + Ftinert
M[:, 0:4] = Ct * Vt * Fss[:, 0:4] / st[:, None]
ss = Fss[:, 4:8].sum(axis=1) + Fs0
M[:, 4:8] = Ct * Vs * Fss[:, 4:8] / ss[:, None]
M = M[1:, :].reshape(-1)
# Relax the transient model to steady state.
for _ in range(n_relax):
sol = spint.solve_ivp(state_space_rhs, [0, 1], M,
args=(Qt, Qs, yCH4, yC2H4, yinert),
rtol=1e-8, atol=1e-8, method='LSODA')
M = sol.y[:, -1]
return M
def dma_mr_cocurrent_step(M: np.ndarray, u: np.ndarray, d: float | None = None,
dt: float = 1.0) -> tuple[np.ndarray, np.ndarray]:
"""
Advance the transient co-current DMA-MR by ``dt`` minutes from holdup state
``M``, with inputs ``u = [tube flowrate, shell flowrate]`` held constant and
disturbance ``d`` (methane feed mole percent). Returns ``(M_next, y)``, the
step contract expected by ``dynamic_operability``.
"""
yCH4, yC2H4, yinert = feed_fractions(d)
Qt, Qs = float(u[0]), float(u[1])
# Integrate the transient balances over the requested time interval.
sol = spint.solve_ivp(state_space_rhs, [0, dt], np.asarray(M, dtype=float),
args=(Qt, Qs, yCH4, yC2H4, yinert),
rtol=1e-6, atol=1e-9, method='LSODA')
M_next = sol.y[:, -1]
return M_next, dma_mr_cocurrent_outputs(M_next)
Then, defining the nominal operating point and the operability set bounds:
# Nominal operating point
Qt_nom = 800.0 # Tube volumetric flowrate [dm³/h]
Qs_nom = 1200.0 # Shell volumetric flowrate [dm³/h]
d_nom = 94.2 # Methane feed mole fraction [%]
# Initial condition: the steady state at the nominal point
x0 = dma_mr_cocurrent_x0([Qt_nom, Qs_nom], d=d_nom)
y0 = dma_mr_cocurrent_outputs(x0)
print(f'Nominal steady state: benzene = {y0[0]:.2f} %, H2 = {y0[1]:.2f} %')
# Available Input Set (AIS): +/- 200 dm³/h around the nominal flowrates
AIS_bounds = np.array([[Qt_nom - 200, Qt_nom + 200],
[Qs_nom - 200, Qs_nom + 200]])
# Desired Output Set (DOS): target product concentrations
DOS_bounds = np.array([[5.0, 9.0], # Benzene mole fraction [%]
[19.0, 24.0]]) # Hydrogen mole fraction [%]
Nominal steady state: benzene = 6.61 %, H2 = 21.96 %
Mapping the dynamic achievable output funnel#
The dynamic Achievable Output Set at time step \(k\), \(AOS^k\), is the set of all output values reachable from the initial state within \(k\) steps by choosing inputs within the AIS [10]. Stacked along the time axis, the \(AOS^k\) slices form the characteristic operability funnel.
Importing opyrability’s dynamic operability module:
from opyrability import dynamic_operability
A single call to dynamic_operability runs the whole analysis: it builds the funnel, evaluates the dOI against the DOS at each time step, and plots the funnel with every time slice colored by its dOI.
To build the funnel, opyrability has to find, at each step \(k\), the outputs the reactor can reach. It chooses how to do this automatically from the model. Here the reactor carries 160 internal states, far too many to track with geometry, so it uses simulation: it runs the model forward from the steady state for many different input choices and collects the outputs reached at each step. This stays practical no matter how many internal states the model has, because we only ever record the two outputs.
To check that the funnel really captures everything reachable, a batch of random input sequences (a Monte Carlo sample) is also simulated and overlaid as dots; every dot should fall inside the funnel.
The figure is interactive: drag to rotate, scroll to zoom, and right-drag to pan, both in Jupyter and on the compiled documentation website.
# Step model with the disturbance held at its nominal value
def step(x, u):
return dma_mr_cocurrent_step(x, u, d=d_nom)
results = dynamic_operability(step, x0, AIS_bounds,
DOS=DOS_bounds,
k_max=8,
y0=y0,
monte_carlo=80,
seed=42,
orientation='vertical',
engine='plotly',
labels=['Benzene mole fraction [%]',
'H2 mole fraction [%]'])
print(f'Method selected automatically: {results["method"]}')
print(f'dOI at each minute: {np.round(results["dOI"], 1)}')
Method selected automatically: n-step simulation
dOI at each minute: [ 0. 0.5 1. 2.2 4.3 7.7 11.4 14.2 17.7]
The funnel starts as a single point (the nominal steady state) and expands as the manipulated inputs move the process over time. The colorbar shows the dOI directly on the funnel: at \(k = 0\) the steady state lies inside the DOS but the achievable set has zero measure, and the dOI grows with each time step as the funnel expands into the DOS, reaching about \(18\%\) of the desired output set after 8 minutes. This result is irrespective of the controller type and quantifies an inherent dynamic property of this modular design.
Assessing robustness to methane feed disturbances#
The shale gas feeding a modular unit varies in composition from well to well. To check whether the desired outputs remain achievable regardless of the feed quality, the funnel is built for two disturbance scenarios bracketing the nominal feed, and their step-by-step intersection is computed: any output in the intersection is achievable no matter which disturbance is realized [10].
Importing opyrability’s scenario evaluation module:
from opyrability import dynamic_operability_scenarios
Each feed composition has its own steady state and its own dynamics, so the analysis has to be repeated for every disturbance value. Rather than set this up by hand, dynamic_operability_scenarios asks for two small helper functions: one that returns the reactor’s starting steady state for a given feed, and one that returns the step model for that feed. It then builds one funnel per scenario and computes (and plots) the region they all share, that is, the outputs reachable no matter which feed occurs:
# For a given methane feed d, build the matching step model and the matching
# starting steady state (each feed composition has its own dynamics and its
# own steady state, so both depend on d).
def make_step(d):
def step(x, u):
return dma_mr_cocurrent_step(x, u, d=d)
return step
def make_initial_state(d):
return dma_mr_cocurrent_x0([Qt_nom, Qs_nom], d=d)
scenario_results = dynamic_operability_scenarios(
make_step, make_initial_state, AIS_bounds,
scenarios={'Low quality well (92% CH4)': 92.0,
'High quality well (96% CH4)': 96.0},
DOS=DOS_bounds,
k_max=6,
orientation='vertical',
engine='plotly',
labels=['Benzene mole fraction [%]',
'H2 mole fraction [%]'])
print('Robust intersection area at each minute:',
np.round(scenario_results['intersection']['volumes'], 3))
Robust intersection area at each minute: [0. 0. 0.002 0.062 0.197 0.478]
Each scenario produces its own funnel (colored outlines), starting from a different steady state because the feed composition shifts the operating point. The blue filled region is their intersection: for the first two minutes the funnels do not overlap and no output can be guaranteed, but from \(k = 3\) onward the intersection grows steadily, meaning adjusting the inlet flowrates can compensate for the feed quality and steer the outputs into a common achievable region.
Conclusions#
Key Results:
The DMA-MR dynamic achievable output funnel was mapped from the nominal steady state with a single call to
dynamic_operability, which automatically chose to map the funnel by simulation, the right approach for this 160-state reactor model.The dOI, displayed directly on the funnel, grows from \(0\%\) to about \(18\%\) of the DOS within 8 minutes of operation.
Under methane feed disturbances of \(92\%\) to \(96\%\), a nonempty robust achievable region appears after 3 minutes, showing that adjusting the inlet flowrates can compensate for the feed-quality disturbance.
Workflow Summary:
Define the dynamic step model and obtain the initial steady state.
Define the AIS and DOS bounds.
Call
dynamic_operabilityto map the funnel and evaluate the dOI.Call
dynamic_operability_scenariosto assess robustness to disturbances.