[Python] How to implement dynamic modeling of a plug flow reactor with catalyst deactivation...

I'm working on a dynamic simulation of a tubular plug flow reactor (PFR) for a gas-phase reaction with time-dependent catalyst deactivation.

The reaction is a first-order irreversible process:

A → B
-r_A = k * C_A
k = k0 * exp(-Ea / (R * T)) Catalyst deactivation follows a power-law model:

dX/dt = -kd * X^n
X(t=0) = 1 I want to simulate this in Python and study the effects of temperature ramps and varying feed composition over time. I'm currently considering scipy.integrate.solve_ivp but also evaluating Pyomo.DAE for more advanced control and optimization possibilities.

My main questions are:

How can I couple the axial spatial domain of the PFR with the time-varying catalyst activity (X(t))?

What's the best approach for implementing this as a system of ODEs or PDEs (e.g. via method of lines)?

Are there efficient discretization schemes or Pyomo/DynPy templates to keep computation tractable?

Bonus: How can I structure the problem to allow for parameter estimation or process optimization later on?

If you’ve tackled similar systems or have experience modeling reactor dynamics with Pyomo or SciPy, I’d appreciate your guidance.

So far, I’ve experimented with scipy.integrate.solve_ivp by discretizing the reactor length using the method of lines. I implemented a system of ODEs where each segment has its own local concentration and deactivation state. The model works for small systems, but becomes stiff and unstable with realistic discretization or temperature profiles.

I also tried using Pyomo.DAE for better control, but I’m still unsure how to correctly couple the axial flow in the PFR with the time evolution of catalyst activity without overcomplicating the structure.

I expected a smoother integration between space and time domains and hoped for a scalable way to simulate more realistic cases (e.g., startup/shutdown, thermal swings, varying inlet conditions).

Continue reading...