A Stability/Bifurcation Framework For Process Design

1. A Stability/Bifurcation Framework For Process Design C. Theodoropoulos1, N. Bozinis2, C. Siettos1, C.C. Pantelides2 and I.G. Kevrekidis11Department of Chemical Engineering,Princeton University, Princeton, NJ 08544 2 Centre for Process System Engineering, Imperial College, London, SW7 2BY, UK

2. bif. quantity parameter Motivation • A large number of existing scientific, large-scale legacy codes • Based on transient (timestepping) schemes. • Enable legacy codes perform tasks such as bifurcation/stability analysis • Efficiently locate multiple steady states and assess the stability of solution branches. • Identify the parametric window of operating conditions for optimal performance • Locate periodic solutions • Autonomous, forced (PSA,RFR) • Appropriate controller design. • RPM: method of choice to build around existing time-stepping codes. • Identifies the low-dimensional unstable subspace of a few “slow” eigenvalues • Stabilizes (and speeds-up) convergence of time-steppers even onto unstable steady-states. • Efficient bifurcation analysis by computing only the few eigenvalues of the small subspace. • Even when Jacobians are not explicitly available (!)

3. Reconstruct solution: u= p+q = PN(p,q)+QF Initial state un Newton iterations Picard iterations Timestepping Legacy Code F(un) Picard iteration NO Subspace Q=I-P Subspace Pof few slow eigenmodes Convergence? YES Final state uf Recursive Projection Method (RPM) • Treats timstepping routine, as a “black-box” • Timestepper evaluates un+1= F(un) F.P.I. • Recursively identifies subspace of slow eigenmodes, P • Substitutes pure Picard iteration with • Newton method inP • Picard iteration in Q = I-P • Reconstructs solution u from sum of the projectors P and Q onto subspace P and its orthogonal complement Q, respectively: • u = PN(p,q) + QF

4. Nonlinear programming solvers Nonlinear algebraic equation solvers gPROMS Steady-state & Dynamic Simulation Steady-state & Dynamic Optimisation Differential algebraic equation solvers Dynamic optimisation solvers gPROMS Model Parameter Estimation Data Reconciliation Maximum likelihood estimation solvers gPROMS:A General Purpose Package

5. Mathematical solution methods in gPROMS • Combined symbolic, structural & numerical techniques • symbolic differentiation for partial derivatives • automatic identification of problem sparsity • structural analysis algorithms • Advanced features: • exploitation of sparsity at all levels • support for mixed analytical/numerical partial derivatives • handling of symmetric/asymmetric discontinuities at all levels • Component-based architecture for numerical solvers • open interface for external solver components • hierarchical solver architectures • mix-and-match • external solvers can be introduced at any level of the hierarchy • well-posedness • DAE index analysis • consistency of DAE IC’s • automatic block triangularisation

6. FitzHugh-Nagumo: An PDE-based Model • Reaction-diffusion model in one dimension • Employed to study issues of pattern formation in reacting systems • e.g. Beloushov-Zhabotinski reaction • u “activator”, v “inhibitor” • Parameters: • no-flux boundary conditions • e, time-scale ratio, continuation parameter • Variation of e produces turning points and Hopf bifurcations

7. Bifurcation Diagrams Around Hopf Around Turning Point <u> e

8. Eigenspectrum Around Hopf

9. Eigenvectors e = 0.02

10. System: y F.P.I p Pseudo – arc length condition (1) (2) Solve (1) & (2) continuation (II) through FORTRAN Arc-length continuation with gPROMS continuation (I) within gPROMS

11. ODEs : DAEs : F.P.I F.P.I R.P.M. through FORTRAN Continuation within gPROMS Getting system Jacobian through an FPI System Jacobian Obtain “correct” Jacobian of leading eigenspectrum Cannot get “correct” Jacobian from augmented system Jacobian of the ODE Stability matrix

12. Tubular Reactor: A DAE system Dimensionless equations: (1) (2) Boundary Conditions: (3) (4) Eqns (1)-(4): system of DAEs. Can also substitute to obtain system of ODEs.

13. Bifurcation/Stability with RPM-gPROMS • Model solved as DAE system • 2 algebraic equations @ each boundary • 101-node FD discretization • 2 unknowns (x1,x2) per node Hopf pt. • State variables: • 99 (x 2) unknowns at inner nodes • Perform RPM-gPROMs at 99-space • to obtain correct Jacobian

14. 1 0.5 0 -1 -0.5 0 0.5 1 1.5 -0.5 -1 Da =0.121738 1 0.5 0 -1 -0.5 0 0.5 1 1.5 -0.5 -1 Da =0.110021 Eigenspectrum

15. Stability Analysis without the Equations SYSTEM AROUND STEADY STATE Leading Spectrum y(k) Matrix-free ARNOLDI + ε q Large-scale eigenvalue calculations (Arnoldi using system Jacobian): R.B. Lechouq & A.G. Salinger, Int. J. Numer. Meth.(2001)

16. t= T/2 to T P(z=0)=Pw z=L Mass balance in ads. bed Darcy’s law z=0 Rate of ads. Rapid Pressure Swing Adsorption1-Bed 2-Step Periodic Adsorption Process t=0 to T/2 Ci(z=0)=PfYf/(RTf) P(z=0)=Pf • Isothermal operation • Modeling Equations (Nilchan & Pantelides) Step 2: Depressurisation Step 1 : Pressurisation

17. Rapid Pressure Swing Adsorption1-Bed 2-Step Periodic Adsorption Process q , c(t=T) • Production of oxygen enriched air • Zeolite 5A adsorbent (300m) • Bed 1m long, 5cm diameter • Short cycle • 1.5s pressurisation, 1.5s depressurisation • T= 3s • Low feed pressure (Pf = 3 bar) • Periodic steady-state operation • reached after several thousand cycles q ,c(t=0) q, c(t=T/2) Must obtain: q , c(t=T) = q , c(t=0)

18. Typical RPSA simulation results(Nilchan and Pantelides, Adsorption, 4, 113-147, 1998) c1(z=0.5) (mol/m3) Time (s)

19. PRM-gPROMS Spatial Profiles (t=T) z z z z

20. Leading Eigenvectors, l=0.99484 c1 q1 q1 c1 q2 c2 c2 q2

21. Conclusions • Can construct a RPM-based computational framework around large-scale timestepping legacy codes to enable them converge to unstable steady states and efficiently perform bifurcation/stability analysis tasks. • gPROMS was employed as a really good simulation tool • communication with wrapper routines through F.P.I. • Both for PDE and DAE-based systems. • Have “brought to light” features of gPROMS for continuation around turning points and information on the Jacobian and/or stability matrix at steady states of systems. • Employed matrix-free Arnoldi algorithms to perform stability analysis of steady state solutions without having either the Jacobian or even the equations! • Used the RPM-based superstructure to speed-up convergence and perform stability analysis of an almost singular periodically-forced system • Have enabled gPROMS to trace autonomous limit cycles • Newton-Picard computational superstructure for autonomous limit cycles.

22. gPROMS • General purpose commercial package for modelling, optimization and control of process systems. • Allows the direct mathematical description of distributed unit operations • Operating procedures can be modelled • Each comprising of a number of steps • In sequence, in parallel, iteratively or conditionally. • Complex processes: combination of distributed and lumped unit operations • Systems of integral, partial differential, ordinary differential and algebraic equations (IPDAEs). • gPROMS solves using method of lines family of numerical methods. • Reduces IPDAES to systems of DAEs. • Time-stepping or pseudo-timestepping. • Jacobians NOT explicitly available. • Cannot perform systematic bifurcation/stability analysis studies.

23. F.P.I F.P.I F.P.I Tracing limit cycles Tracing Limit Cycles continuation (I) within gPROMS R.P.M through FORTRAN continuation (II) through FORTRAN Getting system Jacobian through an FPI tracing limit cycles within gPROMS

24. Tracing limit cycles Tracing Limit Cycles SYSTEM: Periodic Solutions: y(t+T)=y(t) Period T not known beforehand